{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2023-05-27T09:13:08.053227300Z",
     "start_time": "2023-05-27T09:13:04.888934800Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import cv2"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Question2: NumPy Warm-up"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### 1. Choose a few different one-dimensional Gaussian functions (by choosing different mean and variance values), plot them."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [],
   "source": [
    "# plot one-dimensional Gaussian probability distribution function\n",
    "def plot_gaussian_pdf(mu, sigma):\n",
    "    x = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)\n",
    "    y = np.exp(-(x - mu)**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi))\n",
    "    plt.plot(x, y)\n",
    "    plt.show()\n",
    "# choosing different mean and variance values, plot them\n",
    "# plot_gaussian_pdf(0, 1)\n",
    "# plot_gaussian_pdf(0, 2)\n",
    "# plot_gaussian_pdf(0, 3)\n",
    "# plot_gaussian_pdf(1, 1)\n",
    "# plot_gaussian_pdf(1, 2)\n",
    "# plot_gaussian_pdf(1, 3)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-05-27T09:13:08.078881700Z",
     "start_time": "2023-05-27T09:13:08.060233500Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### 2. Verify the above identity for each Gaussian function."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9972920676269242\n",
      "0.9972920676269242\n",
      "0.997292067626924\n",
      "0.9972920676269241\n",
      "0.9972920676269242\n",
      "0.997292067626924\n"
     ]
    }
   ],
   "source": [
    "# Verify the identity ∫ G(x)dx = 1 for each Gaussian function\n",
    "def gaussian_integral(mu, sigma):\n",
    "    x = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)\n",
    "    y = np.exp(-(x - mu)**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi))\n",
    "    return np.trapz(y, x)\n",
    "print(gaussian_integral(0, 1))\n",
    "print(gaussian_integral(0, 2))\n",
    "print(gaussian_integral(0, 3))\n",
    "print(gaussian_integral(1, 1))\n",
    "print(gaussian_integral(1, 2))\n",
    "print(gaussian_integral(1, 3))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-05-27T09:13:08.100447100Z",
     "start_time": "2023-05-27T09:13:08.083576700Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Question3: Numerics and Linear Algebra"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### 1. Use an LU solve to estimate the monomial coefficients $\\vec{c}$. Report the residual L2 norm for both linear systems when N = 8 and N = 16."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2.220446049250313e-16, 3.8459253727671276e-16)\n",
      "(1.2560739669470201e-15, 1.1430445635548515e-15)\n"
     ]
    }
   ],
   "source": [
    "# construct a n*n vandermonde matrix\n",
    "def vandermonde_matrix(n):\n",
    "    x = np.linspace(0, 1, n + 1)[:-1]\n",
    "    return np.vander(x, n, increasing=True)\n",
    "# construct a n*n finite Fourier Series basis matrix\n",
    "# Fi,j−1 =\n",
    "# {\n",
    "#   sin (jπx_i), if 1 ≤ j ≤ N/2\n",
    "#   cos ((j − N/2)πx_i), N/2 + 1 ≤ j ≤ N\n",
    "# }\n",
    "def fourier_basis_matrix(n):\n",
    "    x = np.linspace(0, 1, n + 1)[:-1]\n",
    "    y = np.zeros((n, n))\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            if j + 1 <= n/2:\n",
    "                y[i][j] = np.sin((j + 1) * np.pi * x[i])\n",
    "            else:\n",
    "                y[i][j] = np.cos((j + 1 - n/2) * np.pi * x[i])\n",
    "    return y\n",
    "# print (vandermonde_matrix(8))\n",
    "# print (fourier_basis_matrix(8))\n",
    "# Use an LU solve (scipy.linalg.lu from SciPy package) to estimate the monomial coefficients. Report the residual L2 norm for both linear systems when N = 8 and N = 16.\n",
    "def monomial_coefficients(n):\n",
    "    # f(x) = 1 / (1 + x^2)\n",
    "    x = np.linspace(0, 1, n + 1)[:-1]\n",
    "    y = 1 / (1 + x**2)\n",
    "\n",
    "    # use LU solve to estimate the monomial coefficients\n",
    "    vm = vandermonde_matrix(n)\n",
    "    fm = fourier_basis_matrix(n)\n",
    "    vm_lu = sp.linalg.lu_factor(vm)\n",
    "    fm_lu = sp.linalg.lu_factor(fm)\n",
    "    vm_x = sp.linalg.lu_solve(vm_lu, y)\n",
    "    fm_x = sp.linalg.lu_solve(fm_lu, y)\n",
    "    # print(\"monomial coefficients when N = \" + str(n) + \":\")\n",
    "    # print(\"vandermonde matrix: \" + str(vm_x))\n",
    "    # print(\"fourier basis matrix: \" + str(fm_x))\n",
    "    vm_y = np.dot(vm, vm_x)\n",
    "    fm_y = np.dot(fm, fm_x)\n",
    "    vm_r = np.linalg.norm(vm_y - y)\n",
    "    fm_r = np.linalg.norm(fm_y - y)\n",
    "    return vm_r, fm_r\n",
    "print(monomial_coefficients(8))\n",
    "print(monomial_coefficients(16))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-05-27T09:13:09.703786600Z",
     "start_time": "2023-05-27T09:13:08.100447100Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### 2. Plot $N$ vs. $cond( V )$ and $N$ vs. $cond(F)$ for N = 4, 6, 8, ...32."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[163.42742324166863, 8215.649919892317, 450050.0666487678, 25681357.137469992, 1500717461.2546258, 89069186916.48735, 5344580147535.864, 323245758981552.2, 1.9496127860001332e+16, 1.620597993853218e+18, 1.786397875883056e+19, 6.358601189998436e+18, 8.239713325560794e+19, 7.556656808563421e+18, 4.5887143780777574e+18]\n",
      "[5.02733949212585, 31.635891538074045, 236.58740928194547, 1918.313054823126, 16306.251662479055, 143030.8296855278, 1283150.6173285618, 11708033.664982196, 108251959.76479125, 1011580282.6936167, 9535670909.254139, 90545276003.46182, 865095433988.6373, 8309046218664.75, 80487025496112.38]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGsCAYAAACB/u5dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJYklEQVR4nO3dd3xUVd4G8OdmWnogkEJIgkgvoUkRkI6giwqyK66igPKuDUTEVWERFBERRcUVFttSrAiuKOouiihFpCgdCb2XFFomdWYyc94/kjuZIQnJJPfeuZM8389n1s3kzj0nMjJPTvkdSQghQERERKSAIH93gIiIiGoOBgsiIiJSDIMFERERKYbBgoiIiBTDYEFERESKYbAgIiIixTBYEBERkWIYLIiIiEgxDBZERESkGAYLIiIiUozfgsWGDRtw++23IyEhAZIk4auvvvLp9QUFBRgzZgxSUlJgNBoxbNiwMq9bsGABWrVqhZCQELRo0QIffvhh9TtPREREZfJbsMjNzUX79u2xYMGCKr3e6XQiJCQEEyZMwMCBA8u8ZuHChZgyZQpeeOEF/PHHH5gxYwbGjRuHb775pjpdJyIionJIejiETJIkrFy50mvUwWazYerUqfjss89w5coVtG3bFnPmzEHfvn1LvX7MmDG4cuVKqVGPHj16oGfPnnjttdfczz311FPYunUrfvnlF5V+GiIiotpLt2ssxo8fj82bN2PZsmXYs2cP7rrrLtxyyy04fPhwpe9hs9kQHBzs9VxISAi2bdsGh8OhdJeJiIhqPV0Gi1OnTmHx4sVYsWIFevXqhSZNmuDvf/87brrpJixevLjS9xk8eDA++OADbN++HUII/P777/jggw/gcDhw4cIFFX8CIiKi2sno7w6UZe/evXA6nWjevLnX8zabDfXq1av0faZNm4a0tDTceOONEEIgLi4Oo0ePxquvvoqgIF1mKiIiooCmy2CRk5MDg8GA7du3w2AweH0vPDy80vcJCQnBokWL8O677yI9PR0NGjTAe++9h4iICMTExCjdbSIiolpPl8GiY8eOcDqdyMjIQK9evap9P5PJhMTERADAsmXLcNttt3HEgoiISAV+CxY5OTk4cuSI++vjx49j165diI6ORvPmzTFy5EiMGjUKr7/+Ojp27IjMzEysXbsW7dq1w5AhQwAA+/fvh91ux6VLl5CdnY1du3YBADp06AAAOHToELZt24Zu3brh8uXLeOONN7Bv3z4sXbpU6x+XiIioVvDbdtN169ahX79+pZ4fPXo0lixZAofDgZdeegkffvghzp49i/r16+PGG2/EjBkzkJKSAgC47rrrcPLkyVL3kH+k1NRU3HvvvTh48CBMJhP69euHOXPmoEWLFur+cERERLWULupYEBERUc3AhQZERESkGAYLIiIiUozmizddLhfOnTuHiIgISJKkdfNERERUBUIIZGdnIyEh4Zo7KzUPFufOnUNSUpLWzRIREZECTp8+7S7hUBbNg0VERASAoo5FRkZq3TwRERFVgdVqRVJSkvtzvDyaBwt5+iMyMpLBgoiIKMBUtIyBizeJiIhIMQwWREREpBgGCyIiIlKMLg8hczqdcDgc/u4GBSiTyVTqVFwiItKG7oJFTk4Ozpw5A1Yap6qSJAmJiYkIDw/3d1eIiGodXQULp9OJM2fOIDQ0FDExMSygRT4TQiAzMxNnzpxBs2bNOHJBRKQxXQULh8MBIQRiYmIQEhLi7+5QgIqJicGJEyfgcDgYLIiINKbLxZscqaDq4PuHiMh/dBksiIiIKDAxWBAREZFiGCwCmCRJ+Oqrr/zdDUUtWbIEderU8Xc3iIioihgsqun222/HLbfcUub3Nm7cCEmSsGfPHo17RURE5B8MFtU0duxYrFmzBmfOnCn1vcWLF6Nz585o166dH3pWMbvd7u8uEBEpJs9eiIXrjuL4hVx/d6VW03WwEEIgz17ol0dlC3TddtttiImJwZIlS7yez8nJwYoVKzBs2DDcc889aNiwIUJDQ5GSkoLPPvvM69q+fftiwoQJeOaZZxAdHY34+Hi88MILXtccPnwYvXv3RnBwMFq3bo01a9aU6svp06cxYsQI1KlTB9HR0Rg6dChOnDjh/v6YMWMwbNgwzJo1CwkJCWjRogVOnDgBSZKwfPly9OrVCyEhIejSpQsOHTqE3377DZ07d0Z4eDhuvfVWZGZmuu/lcrnw4osvIjExERaLBR06dMDq1avd35fv++WXX6Jfv34IDQ1F+/btsXnzZq8+L1myBMnJyQgNDcWdd96Jixcvlvq5vv76a3Tq1AnBwcG4/vrrMWPGDBQWFlb0R0NEtcw3u89hzuoDmPvDQX93pVbTVR2Lq+U7nGg9/Xu/tL3/xcEINVf8r8doNGLUqFFYsmQJpk6d6t7quGLFCjidTtx3331YsWIFnn32WURGRuK7777D/fffjyZNmqBr167u+yxduhSTJk3C1q1bsXnzZowZMwY9e/bEzTffDJfLheHDhyMuLg5bt25FVlYWJk6c6NUPh8OBwYMHo3v37ti4cSOMRiNeeukl3HLLLdizZw/MZjMAYO3atYiMjCwVTJ5//nnMmzcPycnJePDBB3HvvfciIiICb731FkJDQzFixAhMnz4dCxcuBAC89dZbeP311/Huu++iY8eOWLRoEe644w788ccfaNasmfu+U6dOxdy5c9GsWTNMnToV99xzD44cOQKj0YitW7di7NixmD17NoYNG4bVq1fj+eef9+rXxo0bMWrUKPzzn/9Er169cPToUTz00EPuPhMRyU5dygMAnC7+J/mHJDSunW21WhEVFYWsrCxERkZ6fa+goADHjx9H48aNERwcjDx7oe6DBQAcOHAArVq1ws8//4y+ffsCAHr37o1GjRrho48+KnX9bbfdhpYtW2Lu3LkAikYsnE4nNm7c6L6ma9eu6N+/P1555RX88MMPGDJkCE6ePImEhAQAwOrVq3Hrrbdi5cqVGDZsGD7++GO89NJLSE1NdYcbu92OOnXq4KuvvsKgQYMwZswYrF69GqdOnXIHjRMnTqBx48b44IMPMHbsWADAsmXLcM8992Dt2rXo378/AOCVV17BkiVLcODAAQBAw4YNMW7cOPzjH//w6nOXLl2wYMGCMu+7f/9+tGnTBqmpqWjZsiXuvfdeZGVl4bvvvnPf469//StWr16NK1euAAAGDhyIAQMGYMqUKe5rPv74YzzzzDM4d+5cmX8eV7+PiKh2+PuK3fhi+xnERliwbepAf3enxrnW57cnXY9YhJgM2P/iYL+1XVktW7ZEjx49sGjRIvTt2xdHjhzBxo0b8eKLL8LpdOLll1/G8uXLcfbsWdjtdthsNoSGhnrd4+p1GA0aNEBGRgYAIDU1FUlJSe5QAQDdu3f3un737t04cuQIIiIivJ4vKCjA0aNH3V+npKS4Q0V57cfFxbmv9XxO7o/VasW5c+fQs2dPr3v07NkTu3fvLve+DRo0AABkZGSgZcuWSE1NxZ133ul1fffu3b2mVHbv3o1NmzZh1qxZ7uecTicKCgqQl5dX6t8jEdVe6dYCAMCFHBsKnS4YDbqe7a+xdB0sJEmq9KiBv40dOxaPP/44FixYgMWLF6NJkybo06cP5syZg7feegvz5s1DSkoKwsLCMHHixFILJ00mk9fXkiTB5XJVuv2cnBzccMMN+OSTT0p9LyYmxv3/w8LCyny9Z/vyiMfVz/nSn2vd19efa8aMGRg+fHip73E0gog8pWUVBQuXADJzbGgQxaMh/CEwPrUDwIgRI/DEE0/g008/xYcffohHH30UkiRh06ZNGDp0KO677z4ARR+qhw4dQuvWrSt971atWuH06dM4f/68+7f+LVu2eF3TqVMnfP7554iNjb3mEJUSIiMjkZCQgE2bNqFPnz7u5zdt2uS1bqQirVq1wtatW72eK+vnOnjwIJo2bVq9ThNRjSePWABFIYPBwj98GidyOp2YNm0aGjdujJCQEDRp0gQzZ87kEecAwsPDcffdd2PKlCk4f/48xowZAwBo1qwZ1qxZg19//RWpqal4+OGHkZ6e7tO9Bw4ciObNm2P06NHYvXs3Nm7ciKlTp3pdM3LkSNSvXx9Dhw7Fxo0bcfz4caxbtw4TJkwocytsdT399NOYM2cOPv/8cxw8eBCTJ0/Grl278MQTT1T6HhMmTMDq1asxd+5cHD58GPPnz/eaBgGA6dOn48MPP8SMGTPwxx9/IDU1FcuWLcNzzz2n9I9ERAEs3+6EtaBkt5hnyCBt+RQs5syZg4ULF2L+/PlITU3FnDlz8Oqrr+Ltt99Wq38BZezYsbh8+TIGDx7sXg/x3HPPoVOnThg8eDD69u2L+Ph4DBs2zKf7BgUFYeXKlcjPz0fXrl3xf//3f15rDgAgNDQUGzZsQHJyMoYPH45WrVph7NixKCgoUGUEY8KECZg0aRKeeuoppKSkYPXq1Vi1apXXjpCK3HjjjXj//ffx1ltvoX379vjhhx9KBYbBgwfj22+/xQ8//IAuXbrgxhtvxJtvvolGjRop/SMRUQC7OkjI0yKkPZ92hdx2222Ii4vDv//9b/dzf/7znxESEoKPP/64UvfwZVcIUVXwfURU+2w5dhF/fa9kKvWRPk0w+daWfuxRzVPZXSE+jVj06NEDa9euxaFDhwAUrdj/5ZdfcOutt5b7GpvNBqvV6vUgIiJS0tUjFpwK8R+fFm9OnjwZVqsVLVu2hMFggNPpxKxZszBy5MhyXzN79mzMmDGj2h0lIiIqjxwkzIYg2J0uToX4kU8jFsuXL8cnn3yCTz/9FDt27MDSpUsxd+5cLF26tNzXTJkyBVlZWe7H6dOnq91pIiIiT2lZNgBAq4SiIXqOWPiPTyMWTz/9NCZPnoy//vWvAIoKKJ08eRKzZ8/G6NGjy3yNxWKBxWKpfk+JiIjKkZ5dFCTaJ0Zh9+krSLMWQAjhrp9D2vFpxCIvLw9BQd4vMRgMVSqcREREpJT04qmPdol1AAB5dieybTys0B98GrG4/fbbMWvWLCQnJ6NNmzbYuXMn3njjDTz44INq9Y+IiKhCacVTH43rhyIy2AhrQSHSswoQGWyq4JWkNJ+Cxdtvv41p06bhscceQ0ZGBhISEvDwww9j+vTpavWPiIjomoQQyLAWrbGIiwxGfFQwrAU5SLMWoFlcRAWvJqX5FCwiIiIwb948zJs3T6XuEBER+eZyngN2Z9GUfGxEMOIig3EoPYc7Q/yER78REVFAk3eA1Aszw2wMQnxksNfzpC0GiwC1bt06SJKEK1eu+LsrRER+Ja+viC0OFPFRwV7Pk7YYLBQwZswYSJJU6nHkyBHV2uzRowfOnz+PqKgo1dogIgoE8o6Q+Mii0gZxxQFDrm1B2uKx6Qq55ZZbsHjxYq/nYmJiVGnL4XDAbDYjPj6+Wvex2+0wm80K9YqIyD/SixduyiMVnArxL32PWAgB2HP98/DxKHiLxYL4+Hivh8FgwPr169G1a1dYLBY0aNAAkydPRmFhyd7q6667rtRi2A4dOuCFF15wfy1JEhYuXIg77rgDYWFhmDVrVplTIb/88gt69eqFkJAQJCUlYcKECcjNzfVqa+bMmRg1ahQiIyPx0EMP+fQzEhHpkXsqJIJTIXqg7xELRx7wcoJ/2v7HOcAcVq1bnD17Fn/6058wZswYfPjhhzhw4AD+9re/ITg42Cs4VMYLL7yAV155BfPmzYPRaMSxY8e8vn/06FHccssteOmll7Bo0SJkZmZi/PjxGD9+vNdIyty5czF9+nQ8//zz1frZiIj0IqM4QMiBQp4KuZBjg8Ppgsmg79+haxp9B4sA8u233yI8PNz99a233ormzZsjKSkJ8+fPhyRJaNmyJc6dO4dnn30W06dPL1XF9FruvfdePPDAA+6vrw4Ws2fPxsiRIzFx4kQAQLNmzfDPf/4Tffr0wcKFC93Hh/fv3x9PPfVUNX5SIiJ9kUcm4orXWNQLM8NkkOBwCmRm25BQJ8Sf3at19B0sTKFFIwf+atsH/fr1w8KFC91fh4WFYdy4cejevbtXrfqePXsiJycHZ86cQXJycqXv37lz52t+f/fu3dizZw8++eQT93NCCLhcLhw/fhytWrWq1H2IiAJNukdxLAAICpIQGxGMs1fykWYtYLDQmL6DhSRVezpCK2FhYWjatKnPrwsKCoK4aj2Hw+Eo8/7XkpOTg4cffhgTJkwo9T3PAFPRfYiIAonD6cLF3OLFm8XBAigavTh7Jd+9Y4S0o+9gEeBatWqF//znP14n7G3atAkRERFITEwEULRz5Pz58+7XWK1WHD9+3Oe2OnXqhP3791cp3BARBaqMbBuEAEwGCXVDS3a5cQGn/3BFi4oee+wxnD59Go8//jgOHDiAr7/+Gs8//zwmTZrkXl/Rv39/fPTRR9i4cSP27t2L0aNHw2Aw+NzWs88+i19//RXjx4/Hrl27cPjwYXz99dcYP3680j8WEZFupHvsCAkKKpl2dteyYLDQHEcsVNSwYUP897//xdNPP4327dsjOjoaY8eOxXPPPee+ZsqUKTh+/Dhuu+02REVFYebMmVUasWjXrh3Wr1+PqVOnolevXhBCoEmTJrj77ruV/JGIiHRFnuqQF27K3LUsOBWiOQYLBSxZsqTc7/Xp0wfbtm0r9/uRkZFYtmyZ13OjR4/2+vrqNRgA0Ldv31LPd+nSBT/88EO5bZ04caLc7xERBaL0q7aayjgV4j+cCiEiooCVdtWOEFmcu/omy3prjcGCiIgCVrq7hsVVIxbu80IKyhz1JfUwWBARUcByT4VcHSyKp0LyHU5YCwpLvY7Uw2BBREQBq+TIdO/Fm8EmA6JCTAB4GJnWdBksOGxF1cH3D1HtkWEtXRxL5jkdQtrRVbCQ6zfY7XY/94QCmfz+qUo9ECIKHDm2QuTYiqY5rl5jAQBx3BniF7rabmo0GhEaGorMzEyYTCafDukiAgCXy4XMzEyEhobCaNTV25uIFCaPRERYjAizlP7vPb54eoS1LLSlq795JUlCgwYNcPz4cZw8edLf3aEAFRQUhOTkZK/D34io5pGPS4+LKj1aAXhMhXDEQlO6ChYAYDab0axZM06HUJWZzWaOdhHVAlcfl341OXBw8aa2dBcsgKLfOIODy06gREREQOnj0q/GEQv/4K91REQUkMorjiVzH0SWxeqbWmKwICKigCQv3ixrqylQUiTrYq4NDqdLs37VdgwWREQUkNKzrz1iER1qhskgQQggI5ujFlphsCAiooBU3pHpsqAgCbERLJKlNQYLIiIKOC6XcI9CXH1kuqd47gzRHIMFEREFnIu5dhS6BCQJqB9e9ogFwLLe/sBgQUREAUcegagfboHJUP5Hmbz+giMW2mGwICKigFPecelXi48qGs1gLQvtMFgQEVHAqajqpiyOUyGaY7AgIqKAU1HVTVk8p0I0x2BBREQBp2SraUVTISVlvYUQqveLGCyIiCgApVVyjYUcPAocLljzC1XvFzFYEBFRAEqv4Mh0WbDJgDqhJgBcwKkVBgsiIgo46ZVcvAnwlFOtMVgQEVFAsRU6cTnPAaDiqRDAc2dIvqr9oiIMFkREFFAyineEmI1BiAoxVXh9PI9P1xSDBRERBRTP4liSJFV4fVwUp0K0xGBBREQBpbI7QmSsZaEtBgsiIgoochXN2Eos3AQ8ynqz+qYmGCyIiCiguI9Lr+SIBQ8i0xaDBRERBZS0SlbdlMkB5GKuHbZCp2r9oiIMFkREFFAqWxxLFh1mhrn4aHV5Rwmph8GCiIgCijtYRFRujYUkSe71GJwOUR+DBRERBQwhRMmukEqOWACsvqklBgsiIgoY1oJCFDhcACq/xgLwqGXBnSGqY7AgIqKAIU9lRIWYEGwyVPp1rGWhHQYLIiIKGOk+FseSNXBX3+TiTbUxWBARUcDwtTiWzF3LglMhqmOwICKigFHVEYt4nheiGQYLIiIKGOnFUxm+7AgBvHeFCCEU7xeVYLAgIqKAIY84xPo4YiFPndgLXbiS51C8X1SCwYKIiAJGRhWnQixGA6LDzAA4HaI2BgsiIgoYciiI83HxZtFruM5CCwwWREQUEAqdLmT6eLKpp3i5rDd3hqiKwYKIiALCxVw7XAIwBEmoF+77iAV3hmiDwYKIiAKCXMMiJtwCQ5Dk8+vjWH1TEwwWREQUEHw9Lv1q7i2nnApRFYMFEREFBF+PS79aHMt6a4LBgoiIAkJVjkv3xIPItMFgQUREAUGuuunLceme5GBxKdcOW6FTsX6RNwYLIiIKCO6pkCoGizqhJpiNRR97GZwOUQ2DBRERBYSqHkAmkyTJ68wQUgeDBRERBQR5N0dVqm7KuDNEfQwWRESke/l2J6wFhQCqvt3U87VcwKkeBgsiItI9OQiEmAyIsBirfB+5rDdHLNTDYEFERLrnudVUknyvuinjQWTqY7AgIiLdS6/Gqaae4jkVojoGCyIi0r3qbjWVcVeI+nwOFmfPnsV9992HevXqISQkBCkpKfj999/V6BsRERGAkuJYVd1qKis5iMwGIUS1+0Wl+bQC5vLly+jZsyf69euH//3vf4iJicHhw4dRt25dtfpHRETkHmGIVShY2AtduJznQHSYudp9I28+BYs5c+YgKSkJixcvdj/XuHFjxTtFRETkKT2resWxZGZjEOqFmXEx1460rAIGCxX4NBWyatUqdO7cGXfddRdiY2PRsWNHvP/++9d8jc1mg9Vq9XoQERH5Ij1b3hVSvcWbgOd0CNdZqMGnYHHs2DEsXLgQzZo1w/fff49HH30UEyZMwNKlS8t9zezZsxEVFeV+JCUlVbvTRERUewgh3GssYiOqN2IBlOwM4QJOdfgULFwuFzp16oSXX34ZHTt2xEMPPYS//e1veOedd8p9zZQpU5CVleV+nD59utqdJiKi2uNKngP2QhcAILaa200Bj1oWLJKlCp+CRYMGDdC6dWuv51q1aoVTp06V+xqLxYLIyEivBxERUWXJIwvRYWZYjIZq3y+eUyGq8ilY9OzZEwcPHvR67tChQ2jUqJGinSIiIpKlKVTDQiav0+BUiDp8ChZPPvkktmzZgpdffhlHjhzBp59+ivfeew/jxo1Tq39ERFTLZbiPS6/+NAjAqRC1+RQsunTpgpUrV+Kzzz5D27ZtMXPmTMybNw8jR45Uq39ERFTLpWUVLdxUbsSCUyFq8vmIuNtuuw233XabGn0hIiIqRd5qqliwKL7P5TwHChxOBJuqv26DSvCsECIi0jW5OJZSwSIqxASLsejjL6N4Gysph8GCiIh0reTIdGXWWEiSxFoWKmKwICIiXZOLYyk1YuF5LwYL5TFYEBGRbjmcLlzMVT5YuGtZcGeI4hgsiIhItzKzbRACMBkkRIcqd2AYp0LUw2BBRES65T4uPSIYQUGSYvflVIh6GCyIiEi3SnaEKLNwU8apEPUwWBARkW6lK1zOW8ay3uphsCAiIt1KU2FHiOf9Mqw2CCEUvXdtx2BBRES65T4nJErZYBEbUXQ/u9OFS7l2Re9d2zFYEBGRbpWcbKrsGguzMQj1w81ebZAyGCyIiEi3lD4y3ZN8Tx5GpiwGCyIi0q0MldZYACU7Q+TTU0kZDBZERKRLObZC5NgKAZSEACXFyUWysvIVv3dtxmBBRES6JE9RRFiMCLMYFb9/PItkqYLBgoiIdEkuXhWr8MJNWUmw4FSIkhgsiIhIl9Kz1dlqKpOnQlh9U1kMFkREpEvyosq4CHWCBadC1MFgQUREuuQu563SiIUcLLLyHShwOFVpozZisCAiIl2Sg4UaO0IAIDLEiGBT0cdgGqdDFMNgQUREuqRW1U2ZJEmcDlEBgwUREemSmsWxZKy+qTwGCyIi0h2XS6h2ZLqneHeRLAYLpTBYEBGR7lzMtaPQJSBJQEyEOlMhAHeGqIHBgoiIdEceragfboHJoN5HFadClMdgQUREupOu8sJNWQNOhSiOwYKIiHQnvXjhplpbTWXu6pss660YBgsiItIdec1DrMrBIt5jKsTlEqq2VVswWBARke7I53eoPWIRE2GBJAGFLoGLuXZV26otGCyIiEh35API1F5jYTIEoX54URtcwKkMBgsiItIdeTGlmjUsZO4tp1zAqQgGCyIi0p2M7OLFmyodQOYpjrUsFMVgQUREumIrdOJS8XoHtY5M9xQfxakQJTFYEBGRrshnhJiNQagTalK9PU6FKIvBgoiIdMWzOJYkSaq3x6kQZTFYEBGRrsgf8GpvNZXFR7Gst5IYLIiISFfSNTgu3ROnQpTFYEFERLqixXHpnuSy3taCQuTbnZq0WZMxWBARka6kaVR1UxZhMSLUbChqm9Mh1cZgQUREupLuPidE3aqbMkmSOB2iIAYLIiLSlXSNF28CJdMuXMBZfQwWRESkG0KIkiPTNai6KZPb4lRI9TFYEBGRblgLCpHvKFpAqdXiTc+2OBVSfQwWRESkG/JURFSICcEmg2btxkeyrLdSGCyIiEg3PKtuaolTIcphsCAiIt3Q8rh0T+7Fm5wKqTYGCyIi0g33cekaBwt5xCIj2waXS2jadk3DYEFERLrhrxGLmHALgiSg0CVwIdemads1DYMFERHphnuNhYZbTQHAaAhC/fDiBZxZDBbVwWBBRES64Q4WEdou3gS4gFMpDBZERKQb7iPTNR6xADxqWTBYVAuDBRER6YLTJZCZre2R6Z7iuTNEEQwWRESkCxdybHAJwBAkudc7aIlTIcpgsCAiIl2Q11fEhFtgCJI0b58HkSmDwYKIiHShZKup9qMVAHh0ukIYLIiISBdKynlrv74CAOKjigINp0Kqh8GCiIh0QT4u3V/BQm43u6AQefZCv/ShJmCwICIiXfDnVlMAiAg2IcxcdKIqp0OqjsGCiIh0wd9TIUBJxU9Oh1QdgwUREemCv45M9xTPnSHVxmBBRES6IE8/aH2yqaeSnSE8L6SqGCyIiMjv8u1OWAuKFkzG6mAqhCMWVcdgQUREfid/kIeYDIgMNvqtH6xlUX0MFkRE5HfpHjtCJEn7qpsyHkRWfQwWRETkd/IHeawfjkv3FM+pkGpjsCAiIr9L93MNC5k8FZKRbYPTJfzal0DFYEFERH7n76qbsvrhZgRJRUe4X8zhzpCqYLAgIiK/S9NBcSwAMBqCEBPBM0Oqg8GCiIj8LsPq/xoWMu4MqR4GCyIi8rs0HVTdlMWx+ma1VCtYvPLKK5AkCRMnTlSoO0REVNsIIXSzxgIoWUDKqZCqqXKw+O233/Duu++iXbt2SvaHiIhqmSt5DtgLXQCAWB2NWJznVEiVVClY5OTkYOTIkXj//fdRt25dpftERES1iDwyEB1mhsVo8HNveBBZdVUpWIwbNw5DhgzBwIEDK7zWZrPBarV6PYiIiGTpOimOJXNPhXDEokp8Lsi+bNky7NixA7/99lulrp89ezZmzJjhc8eIiKh20EtxLFnJ4k3WsagKn0YsTp8+jSeeeAKffPIJgoMr9waYMmUKsrKy3I/Tp09XqaNERFQzyUeU62GrKVAScHJshcixFfq5N4HHpxGL7du3IyMjA506dXI/53Q6sWHDBsyfPx82mw0Gg/f8mMVigcWij+EtIiLSn/Ts4qkQnQSLcIsR4RYjcmyFSMsqQNPYcH93KaD4FCwGDBiAvXv3ej33wAMPoGXLlnj22WdLhQoiIqKKpGfppziWLC7SgpzMQqRbGSx85VOwiIiIQNu2bb2eCwsLQ7169Uo9T0REVBnyiIUeimPJ4qOCcTQzlws4q4CVN4mIyK/kNRZ6KI4lk/vCIlm+83lXyNXWrVunQDeIiKg2cjhduJhbvHhTJ7tCANayqA6OWBARkd9kZtsgBGAySIgONfu7O26sZVF1DBZEROQ3ae7iWMEICpL83JsSPIis6hgsiIjIb+Tj0vVwRoinBjyIrMoYLIiIyG/SdLjVFCjpT2a2DYVOl597E1gYLIiIyG/Ss/W3IwQA6oVbYAiS4BLAhRy7v7sTUBgsiIjIb+TiWHoLFoYgyX0oGqdDfMNgQUREfpPmPoBMX2ssAI9aFtwZ4hMGCyIi8ht510VchL5GLADWsqgqBgsiIvIb+WjyOB0Vx5LFc2dIlTBYEBGRX3geS663NRaARy0LToX4hMGCiIj8Qp5ikI8p1xt53QdHLHzDYEFERH5RsiNEfws3AR5EVlUMFkRE5Bclx6XrbxoE8Fi8yakQnzBYEBGRX8jHpeut6qZMXryZa3ciu8Dh594EDgYLIiLyi3T3OSH6DBahZiMigovWfnDLaeUxWBARkV/IH9bxOl1jAZSMpsijK1QxBgsiIvKLkqqb+hyxAFjLoioYLIiIyC8yiotj6XUqBPCoZcFgUWkMFkREpDmXS3hMheg3WMTzvBCfMVgQEZHmLuXZUegSkCQgJkK/ayziOBXiMwYLIiLSnDwCUC/MApNBvx9FPIjMd/r90yQiohorXcfHpXviVIjvGCyIiEhz7lNNdXhcuqe44uBzIceGQqfLz70JDAwWRESkOXnNgh6PS/dUP8wCY5AElwAyc1jLojIYLIiISHMZcrDQ+YhFUJCE2OLFpZwOqRwGCyIi0lxagKyxAEpGVbiAs3IYLIiISHNp7iPT9T1iAXABp68YLIiISHMZ2cWLNwMgWMh9TLNyjUVlMFgQEZGmbIVOXMq1A9B31U1ZPKdCfMJgQUREmpLPCDEbg1An1OTn3lSMUyG+YbAgIiJNyb/5x0VaIEmSn3tTMR5E5hsGCyIi0pRcHCsQpkEA76PThRB+7o3+MVgQEZGm5K2mej4u3ZMcgPLsTmTbCv3cG/1jsCAiIk0FwnHpnkLMBkQGGwEA6VxnUSEGCyIi0pTnGotAEc/j0yuNwYKIiDQVSMWxZHHcGVJpDBZERKSpQCqOJYvnzpBKY7AgIiLNCCHcv/UHyhoLgFMhvmCwICIizVgLCpHvcAIIrBGLkqkQlvWuCIMFERFpRj4uPTLYiBCzwc+9qTxOhVQegwUREWmm5Lj0wBmtADgV4gsGCyIi0oxcdTOQpkGAkv5eyLHB4XT5uTf6xmBBRESaKalhEVjBol6YGSaDBCGAzGyus7gWBgsiItJMIO4IAYCgIAmxEZwOqQwGCyIi0kwgVt2UyX1mWe9rY7AgIiLNBOpUCMAFnJXFYEFERJoJ1MWbgEctCwaLa2KwICIiTThdApk5RcEi0LabAh61LDgVck0MFkREpIkLOTY4XQJBElA/PPDWWHAqpHIYLIiISBPy+oqYCAsMQZKfe+M7nnBaOQwWRESkiUDdaiqL91hjIYTwc2/0i8GCiIg0kV5cWCo2UINF8VRIgcMFa36hn3ujXwwWRESkifQAH7EINhkQFWICwHUW18JgQUREmgjUA8g8xXPLaYUYLIiISBPy4s3YiMDbESKLi+KW04owWBARkSbSa8SIRVEo4ohF+RgsiIhIE4FcdVPGqZCKMVgQEZHqChxOZOU7AAR2sOBUSMUYLIiISHVyDYtgUxAig41+7k3VccSiYgwWRESkOvf6ishgSFLgVd2UyaMt6QwW5WKwICIi1aUF8HHpnhoUT4VcyLHDXujyc2/0icGCiIhUl1EDFm4CQHSYGWZD0UdnRjZHLcrCYEFERKqrCcWxAECSJMQWbznldEjZGCyIiEh1NaE4lsy9gDPL5uee6BODBRERqa4mFMeSyVtOuTOkbAwWRESkujSPXSGBLp47Q66JwYKIiFQlhKgRVTdlJVMhDBZlYbAgIiJVXclzuLdmygsfAxmnQq6NwYKIiFSVXrwts26oCRajwc+9qT5OhVwbgwUREalKnjKoCdMggPdUiBDCz73RHwYLIiJSVU3aEQKUTOfYCl3ug9WohE/BYvbs2ejSpQsiIiIQGxuLYcOG4eDBg2r1jYiIagD3ws2ImhEsgk0G1A01AeA6i7L4FCzWr1+PcePGYcuWLVizZg0cDgcGDRqE3NxctfpHREQBzn1OSA0ZsQBKpnW4M6Q0n86uXb16tdfXS5YsQWxsLLZv347evXsr2jEiIqoZMtwHkAX+jhBZfFQwDqRlcwFnGXwKFlfLysoCAERHR5d7jc1mg81WUvbUarVWp0kiIgowNak4loxlvctX5cWbLpcLEydORM+ePdG2bdtyr5s9ezaioqLcj6SkpKo2SUREAUj+8K0pu0IAj6kQjliUUuVgMW7cOOzbtw/Lli275nVTpkxBVlaW+3H69OmqNklERAHG4XThYm7NCxbyDhdOhZRWpamQ8ePH49tvv8WGDRuQmJh4zWstFgsslpozr0ZERJWXmW2DEIAxSEK9MLO/u6MYlvUun0/BQgiBxx9/HCtXrsS6devQuHFjtfpFREQ1gOdx6UFBkp97o5w4Vt8sl0/BYty4cfj000/x9ddfIyIiAmlpaQCAqKgohISEqNJBIiIKXOk1cKspUDIVcjHXDluhs0aUKleKT2ssFi5ciKysLPTt2xcNGjRwPz7//HO1+kdERAHMXc67hhTHktUNNcFsLPoIzbByZ4gnn6dCiIiIKis9u+hDt6aU85ZJkoS4SAtOX8pHurUASdGh/u6SbvCsECKiGuLM5TxsOJSpq18C02vYAWSe4rnltEwMFkRENUCh04X7/70NoxZtw7Lf9LOtXz4yvSZV3ZSxrHfZGCyIiGqAVbvP4fiFonObXvp2P05fyvNzj4rIH7o1qeqmLJ47Q8rEYEFEFOCcLoEFPx8BAISZDci1O/HU8t1wuvw/JSKfbBpbE4NFlDwVwsWbnhgsiIgC3P/2ncfRzFxEhZjwxaM9EGY2YNuJS1j0y3G/9ivHVogcWyGAmrd4E/CoZcGpEC8MFkREAczlEpj/U9FoxQM9r0OrBpGYfntrAMBr3x/EwbRsv/VNniIItxgRbqnWmZe6VDJiwWDhicGCiCiA/ZiajgNp2Qi3GPFAj6JqyCM6J6F/y1jYnS5MWr4L9kKXX/rmrrpZAxduAt67QvS0E8ffGCyIiAKUEALzi9dWjOreCFGhJgBFNRZe+XMK6oaa8Mc5K+b/dNgv/Uuvgcele5IDk73QhSt5Dj/3Rj8YLIiIAtT6Q5nYcyYLISYDxt7kfXZTbEQwZt2ZAgBYsO4odp66rHn/auJx6Z4sRgOiiw9W43RICQYLIqIAJITA28VrK0Z2S0a98NLTDX9KaYBhHRLgdAk8tXw38u1OTfvoPiekhgYLwKOWBYOFG4MFEVEA2nLsErafvAyzMQgP9b6+3Otm3NEW8ZHBOHYhF3NWH9Cwh55TITVzjQVQ8rNxZ0gJBgsiogD0dvG6ib92SbpmjYioUBNe/Us7AMCSX09g05ELmvQPqB0jFtwZUhqDBRFRgNl+8hJ+PXoRJoOEh/s0qfD63s1jcP+NjQAAf1+xG1n52iw0lItj1bQj0z3FsfpmKQwWREQBRl5b8edOiWhYJ6RSr5nyp5a4rl4ozmcVYMY3f6jZPQBF9TUysmvBiAXPCymFwYKIKIDsPZOFdQczYQiS8GjfikcrZKFmI14f0R5BEvDljrNYvS9NxV4Cl/LscDgFJAmIjai5ayziWNa7FAYLIqIAIq+tGNo+AY3qhfn02hsaReOR4qmTqSv34kKOeh+G8m/w9cIsMBlq7kcNDyIrreb+aRMR1TAH0qz4YX86JAl4rF/TKt1j4sDmaNUgEhdz7Zjy5V7VKkZm1ODj0j3JweJSrh22Qm238+oVgwURUYCQzwT5U0oDNI0Nr9I9zMYgvDGiPUwGCWv2p+OL7WeU7KKbXByrplbdlNUJNcFsLPoozeB0CAAGCyKigHA0Mwff7T0PABhfxdEKWasGkZh0cwsAwIxv9uPM5bxq9+9qJeeE1OxgIUmS15khxGBBRBQQFvx8BEIAN7eOQ6sGkdW+30O9r8cNjeoix1aIp1fsgcul7JRITT8nxBN3hnhjsCAi0rlTF/Pw9a5zAIDH+1dvtEJmCJLwxoj2CDUbsPnYRSz59YQi95WlWWvHGgugZGcIF3AWYbAgItK5heuPwOkS6NM8Bu0S6yh230b1wvCPP7UCAMxZfQBHMrIVu3dtKI4lk8t6c8SiCIMFEZGOnbuS715gqdRohaeR3ZLRu3kMbIUuTFq+Gw6nS5H71qapEB5E5o3BgohIx95dfxQOp0D36+uh83XRit9fkiS8+ud2iAoxYc+ZLPzr56PVvqet0IlLuXYANbvqpiyeUyFeGCyIiHQqI7sAn/12GoA6oxWy+KhgzBzWFkBRAa49Z65U637ytkuzIQh1Q03V7Z7ucVeINwYLIiKden/DMdgLXbihUV10b1JP1bbuaJ+AIe0aoNAlMGn5bhQ4ql7sqWSrqQWSJCnVRd1yH0SWZVOt4FggYbAgItKhS7l2fLzlFABgfP+mmnxAvzS0LWIiLDiSkYPXvj9Y5fvICzdrw/oKoCRY2J0u9xRQbcZgQUSkQ//+5RjyHU6kNIxC3+YxmrRZN8yMV//cDgCwaNNxbD56sUr3cW81rQU7QoCiaqb1wswAOB0CMFgQEelOVp4DS389CUC70QpZv5axuKdrMoQA/r5iN7ILHD7fI0MOFhG1I1gAHtMhDBYMFkREerPk1xPIsRWiRVwEbm4Vp3n7U4e0QlJ0CM5eycfMb/f7/Hr5t/b4qJpfHEsm7wyRz0ipzRgsiIh0JMdWiEWbjgMoGq0ICtJ+8WO4xYjX7+oASQKW/34GP+5P9+n1cqGo2rDVVMZaFiUYLIiIdOSjzSeRle/A9TFh+FNKA7/1o2vjaDzU63oAwOQv9+BiTuV/E8/ILq66WYuCRbx7ZwiDBYMFEZFO5Nud+GDjMQDAuL5NYfDDaIWnJ29ujhZxEbiQY8fUlfsqtZVSCOEesagtu0KAkmkfjlgwWBAR6can207hYq4dSdEhuKNDgr+7g2CTAa+PaA9jkITVf6Thq11nK3xNtq0Q+cU1MGrTiAUXb5ZgsCAi0oEChxPvbSgqp/1Y36YwGfTx13PbhlGYOLAZAGD613/gfFb+Na+XpwIig40IMRtU759euBdvMlgwWBAR6cGK7WeQbrWhQVQwhndq6O/ueHmkTxN0SKqD7IJCPL1iD1yu8qdESo5Lrz2jFUDJtM+VPEe1qpbWBAwWRER+5nC68M66otGKR/o0gcWor9/0jYYgvDGiPYJNQfjlyAV8vPVkude6q27WkuJYsqgQE4JNRR+ptX06hMGCiMjPVu44i7NX8lE/3IK7uyT5uztluj4mHFNubQUAePm/qTiWmVPmdem1dMRCkqSSw8hq+c4QBgsiIj8qdLrwr3VHAAAP974ewSZ9jVZ4uv/GRripaX0UOFx4asVuFDpdpa4pCRa1pziWjLUsijBYEBH50bd7zuPExTzUDTXh3m7J/u7ONQUFSXj1L+0QEWzEzlNX8O6GY6WuqY1bTWXy9A+nQoiIyC9cLoH5PxeNVvxfr+sRZjH6uUcVS6gTgheHtgEAvLnmEPadzfL6fsmR6bUwWLinQmp3WW8GCyIiP1n9RxqOZOQgMtiI+7s38nd3Km1Yh4a4pU08Cl0CTy3fDVthyS6I2nZkuifWsijCYEFE5AdCCLz9U9FoxZiejREZbPJzjypPkiTMurMt6oebcTA9G2+sOQQAcLoEMnNq564QoORnPnslv1JVSmsqBgsiIj9Ym5qB1PNWhJkNeLDndf7ujs/qhVswe3g7AMB7G47htxOXcDHHBqdLIEgC6oWZ/dxD7TUoDha7Tl/B4Hkb8N6Go8jIrn2jFwwWREQaE0Lg7eK1Ffd3vw51QgPzQ/jm1nEY0TkRQgCTlu/C0cxcAEBMhAVGnVQO1VJKwyj8tUsSLMYgHErPwcv/PYDus3/Cg0t+w3/3nveaMqrJJKHxeI3VakVUVBSysrIQGRmpZdNERLqw4VAmRi3aVlRw6tn+qB8euFszswscuGXeRpy9ko/r64fh2IVctEuMwqrxN/m7a35jLXDguz3nseL309hx6or7+TqhJgxtn4C/3JCEtg0jIUn+PWTOV5X9/Nb/EmQiohpmfvHainu7NgroUAEAEcEmzL2rPe55fwuOXSgasahtxbGuFhlswj1dk3FP12QczczBf7afwZc7ziLNWoClm09i6eaTaBEXgb/ckIhhHRsiJiKw3wNXq31jVUREfrTl2EVsO3EJZkMQHup9vb+7o4juTeph7E2N3V/XxuJY5WkSE45nbmmJTZP748MHu+L29gkwG4NwMD0bs/6bihtnr8X/Lf0Nq/edh72wdMGxQMQRCyIiDcmjFSO6JNaonRNPD26B9YcycSQjB8nRof7uju4YgiT0bh6D3s1jkJXvwLd7zuGL7Wew89QV/JiagR9TM1A31IShHRriLzckok1C4E2VyLjGgohIIztOXcbwf/0KY5CEdU/3RWLdmvUBfPpSHlbuPIvR3a9DVGjgbJ/1pyMZOfjPjjP4cscZdw0QAGgZXzJVopfpssp+fjNYEBFp5MElv+GnAxkY0TkRr/6lvb+7QzridAlsPJyJL7afwQ/7093TIsYgCX1bxOKuzono1yIWZqP/VjBw8SYRkY7sO5uFnw5kIEgCHu3b1N/dIZ0xFAeIvi1ikZXnwDd7zmHF9jPYffoKfkxNx4+p6YgOM2Noh4TiqZIof3e5XByxICLSwCMfbcfqP9IwtEMC3vprR393hwLE4fRsfLGjaFdJZnbJVEmrBpFFUyUdElBPo6kSToUQEenEofRsDHpzAwBgzZO90Swuws89okBT6HRh45EL+OL3M1izPx12Z8lUSf+WsfjLDYno1zIWJhULk3EqhIhIJ+SdILe2jWeooCoxGoLQr0Us+rWIxZU8O77ZXbSrZPeZLPywPx0/7E9HvTAzhnZoiLs6J6JVA//94s4RCyIiFR3LzMHAN9bDJYDvJtyk67lxCjyH0rPxRXEBrgs5JVMlXzzSHZ2vi1a0LY5YEBHpwL/WHYVLAANaxjJUkOKax0XgH39qhWcGt8CG4l0l+89Z0Sm5rt/6xGBBRKQSua4DAIzvz50gpB6jIQj9W8ahf8s4OJwuBAX5r7gWS3oTEalk4fqjcLoEejWrj45+/A2Sahc1F3BWBoMFEZEKzmfl44vfzwAAHu/fzM+9IdIOgwURkQreXX8MdqcLXRtHo2tjZRfREekZgwURkcIys234bNspAMAEjlZQLcNgQUSksA82HoOt0IUOSXXQs2k9f3eHSFMMFkRECrqca8dHW04CACYMaBqwR18TVRWDBRGRghZtOo48uxNtEiLRr0Wsv7tDpDkGCyIihWTlO7Bk0wkAwOP9OVpBtRODBRGRQpb+egLZtkI0jwvHoNbx/u4OkV+w8iYRURVczLFh79ks7D2ThT3F/0yzFgAAxvVr6tfKh0T+xGBBRFSBK3l27D2bhT1nigLE3rNZOHslv9R1kgT8KaUBbmuX4IdeEukDgwURkQdrgQP7rhqJOHUpr8xrr48JQ7uGUUhJrIOUhlFokxCJMAv/WqXarUr/BSxYsACvvfYa0tLS0L59e7z99tvo2rWr0n0jIlJVjq0Qf5zNKhmNOJuF4xdyy7y2Ub1QpDSMQrvEKKQ0rIM2DSMRGWzSuMdE+udzsPj8888xadIkvPPOO+jWrRvmzZuHwYMH4+DBg4iN5dYqItKnPHsh9p+zeq2LOJqZAyFKX5tYN8QdINolRqFtQhSiQhkiiCpDEqKs/6zK161bN3Tp0gXz588HALhcLiQlJeHxxx/H5MmTK3y91WpFVFQUsrKyEBkZWbVeExFdQ4HDidTzVq91EYczsuEq42+7BlHBJSMRxVMa0WFm7TtNpHOV/fz2acTCbrdj+/btmDJlivu5oKAgDBw4EJs3by7zNTabDTabzatjatjywSTAnq3KvYkChk+/JqhHeHTE81eXa3XP+zrhvrjM1wjP54X7f4UoqiVxOc/uvl/r4gcMQIjZgPrhluKHGfXCLQg1GYoutAE4WvwgCnT9/gEE++eXd5+CxYULF+B0OhEXF+f1fFxcHA4cOFDma2bPno0ZM2ZUvYeV1OTMl4jBZdXbIaIAYCjneRcAa/GDqCa76cnACBZVMWXKFEyaNMn9tdVqRVJSkuLtHGk8EkfsOYrfl4h8I5X1/z0qUHp9/6pSD9I1rr36esnzfz2ejwg2IjYiGOEWY6n7E9Ua5lC/Ne1TsKhfvz4MBgPS09O9nk9PT0d8fNlV5iwWCywWS9V7WEndR89SvQ0iIiK6Np9KepvNZtxwww1Yu3at+zmXy4W1a9eie/fuineOiIiIAovPUyGTJk3C6NGj0blzZ3Tt2hXz5s1Dbm4uHnjgATX6R0RERAHE52Bx9913IzMzE9OnT0daWho6dOiA1atXl1rQSURERLWPz3Usqot1LIiIiAJPZT+/eWw6ERERKYbBgoiIiBTDYEFERESKYbAgIiIixTBYEBERkWIYLIiIiEgxDBZERESkGAYLIiIiUgyDBRERESlG9WPTryYX+rRarVo3TURERFUkf25XVLBb82CRnZ0NAEhKStK6aSIiIqqm7OxsREVFlft9zc8KcblcOHfuHCIiIiBJkmL3tVqtSEpKwunTp/1yBom/29dDH9h+7W5fD31g+7W7fT30oSa3L4RAdnY2EhISEBRU/koKzUcsgoKCkJiYqNr9IyMj/Xq4mb/b10Mf2H7tbl8PfWD7tbt9PfShprZ/rZEKGRdvEhERkWIYLIiIiEgxNSZYWCwWPP/887BYLLWyfT30ge3X7vb10Ae2X7vb10Mfanv7gB8WbxIREVHNVWNGLIiIiMj/GCyIiIhIMQwWREREpBgGCyIiIlJMjQoWr7zyCiRJwsSJEzVt9+zZs7jvvvtQr149hISEICUlBb///rsmbTudTkybNg2NGzdGSEgImjRpgpkzZ1ZYy706NmzYgNtvvx0JCQmQJAlfffWV1/eFEJg+fToaNGiAkJAQDBw4EIcPH9akfYfDgWeffRYpKSkICwtDQkICRo0ahXPnzmnS/tUeeeQRSJKEefPmadp+amoq7rjjDkRFRSEsLAxdunTBqVOnNGk/JycH48ePR2JiIkJCQtC6dWu88847irQNALNnz0aXLl0QERGB2NhYDBs2DAcPHvS6pqCgAOPGjUO9evUQHh6OP//5z0hPT9ek/UuXLuHxxx9HixYtEBISguTkZEyYMAFZWVmatO9JCIFbb721wvepWn3YvHkz+vfvj7CwMERGRqJ3797Iz8/XpP20tDTcf//9iI+PR1hYGDp16oT//Oc/1W4bABYuXIh27dq5i1B1794d//vf/9zfV/P9V1H7ar//KqPGBIvffvsN7777Ltq1a6dpu5cvX0bPnj1hMpnwv//9D/v378frr7+OunXratL+nDlzsHDhQsyfPx+pqamYM2cOXn31Vbz99tuqtZmbm4v27dtjwYIFZX7/1VdfxT//+U+888472Lp1K8LCwjB48GAUFBSo3n5eXh527NiBadOmYceOHfjyyy9x8OBB3HHHHYq0XVH7nlauXIktW7YgISFBsbYr0/7Ro0dx0003oWXLlli3bh327NmDadOmITg4WJP2J02ahNWrV+Pjjz9GamoqJk6ciPHjx2PVqlWKtL9+/XqMGzcOW7ZswZo1a+BwODBo0CDk5ua6r3nyySfxzTffYMWKFVi/fj3OnTuH4cOHa9L+uXPncO7cOcydOxf79u3DkiVLsHr1aowdO1aT9j3NmzdP0aMTfOnD5s2bccstt2DQoEHYtm0bfvvtN4wfP/6apaCVbH/UqFE4ePAgVq1ahb1792L48OEYMWIEdu7cWe32ExMT8corr2D79u34/fff0b9/fwwdOhR//PEHAHXffxW1r/b7r1JEDZCdnS2aNWsm1qxZI/r06SOeeOIJzdp+9tlnxU033aRZe1cbMmSIePDBB72eGz58uBg5cqQm7QMQK1eudH/tcrlEfHy8eO2119zPXblyRVgsFvHZZ5+p3n5Ztm3bJgCIkydPatb+mTNnRMOGDcW+fftEo0aNxJtvvql42+W1f/fdd4v77rtPlfYq036bNm3Eiy++6PVcp06dxNSpU1XpQ0ZGhgAg1q9fL4Qoer+ZTCaxYsUK9zWpqakCgNi8ebPq7Zdl+fLlwmw2C4fDoVn7O3fuFA0bNhTnz5+v1H8nSvehW7du4rnnnlOtzYraDwsLEx9++KHXddHR0eL9999XpQ9169YVH3zwgebvv6vbL4ua77+y1IgRi3HjxmHIkCEYOHCg5m2vWrUKnTt3xl133YXY2Fh07NgR77//vmbt9+jRA2vXrsWhQ4cAALt378Yvv/yCW2+9VbM+eDp+/DjS0tK8/iyioqLQrVs3bN682S99ysrKgiRJqFOnjibtuVwu3H///Xj66afRpk0bTdr0bPu7775D8+bNMXjwYMTGxqJbt26KDoNXpEePHli1ahXOnj0LIQR+/vlnHDp0CIMGDVKlPXmINzo6GgCwfft2OBwOr/dgy5YtkZycrMp78Or2y7smMjISRqPyxzOV1X5eXh7uvfdeLFiwAPHx8Yq3WVEfMjIysHXrVsTGxqJHjx6Ii4tDnz598Msvv2jSPlD0Pvz8889x6dIluFwuLFu2DAUFBejbt6+ibTudTixbtgy5ubno3r275u+/q9svi5rvvzJpEl9U9Nlnn4m2bduK/Px8IYTQfMTCYrEIi8UipkyZInbs2CHeffddERwcLJYsWaJJ+06nUzz77LNCkiRhNBqFJEni5Zdf1qRtIUr/xrpp0yYBQJw7d87rurvuukuMGDFC9favlp+fLzp16iTuvfdexdsur/2XX35Z3HzzzcLlcgkhhKYjFvJvp6GhoeKNN94QO3fuFLNnzxaSJIl169ap3r4QQhQUFIhRo0YJAMJoNAqz2SyWLl2qeNtCFL3/hwwZInr27Ol+7pNPPhFms7nUtV26dBHPPPOM6u1fLTMzUyQnJ4t//OMfirZ9rfYfeughMXbsWPfXFf13onQfNm/eLACI6OhosWjRIrFjxw4xceJEYTabxaFDh1RvXwghLl++LAYNGuR+H0ZGRorvv/9esXb37NkjwsLChMFgEFFRUeK7774TQmj3/iuv/aup+f4rj+anmyrp9OnTeOKJJ7BmzRrF5o995XK50LlzZ7z88ssAgI4dO2Lfvn145513MHr0aNXbX758OT755BN8+umnaNOmDXbt2oWJEyciISFBk/b1zOFwYMSIERBCYOHChZq0uX37drz11lvYsWOHKnPbFXG5XACAoUOH4sknnwQAdOjQAb/++iveeecd9OnTR/U+vP3229iyZQtWrVqFRo0aYcOGDRg3bhwSEhIUH1UcN24c9u3bp9pvwtVt32q1YsiQIWjdujVeeOEFTdpftWoVfvrpJ0XWElS1D/L78OGHH8YDDzwAoOjvxrVr12LRokWYPXu2qu0DwLRp03DlyhX8+OOPqF+/Pr766iuMGDECGzduREpKSrXbbdGiBXbt2oWsrCx88cUXGD16NNavX1/t+1a3/datW7uvUfv9Vy7NIowKVq5cKQAIg8HgfgAQkiQJg8EgCgsLVe9DcnKy128GQgjxr3/9SyQkJKjethBCJCYmivnz53s9N3PmTNGiRQtN2sdVvwkdPXpUABA7d+70uq53795iwoQJqrcvs9vtYtiwYaJdu3biwoULirdbXvtvvvmm+/3n+Z4MCgoSjRo1Ur19m80mjEajmDlzptd1zzzzjOjRo4fq7efl5QmTySS+/fZbr+vGjh0rBg8erGjb48aNE4mJieLYsWNez69du1YAEJcvX/Z6Pjk5Wbzxxhuqty+zWq2ie/fuYsCAAe4RVSWV1/4TTzxR7nuwT58+mvTh2LFjAoD46KOPvJ4fMWKEoqOH5bV/5MgRAUDs27fP6/kBAwaIhx9+WLH2r773Qw89pNn7r7z2ZWq//64loNdYDBgwAHv37sWuXbvcj86dO2PkyJHYtWsXDAaD6n3o2bNnqW1Ohw4dQqNGjVRvGyiaS716lbXBYHD/xqC1xo0bIz4+HmvXrnU/Z7VasXXr1nLn/5Qmj1QcPnwYP/74I+rVq6dJuwBw//33Y8+ePV7vyYSEBDz99NP4/vvvVW/fbDajS5cufntPOhwOOBwOVd+TQgiMHz8eK1euxE8//YTGjRt7ff+GG26AyWTyeg8ePHgQp06dUuQ9WFH7QNF7ftCgQTCbzVi1apWiI6oVtT958uRS70EAePPNN7F48WJN+nDdddchISFBtfdhRe3n5eUBgKZ/N7pcLthsNtXffxW1D6j7/qsUTWOMBrReY7Ft2zZhNBrFrFmzxOHDh8Unn3wiQkNDxccff6xJ+6NHjxYNGzYU3377rTh+/Lj48ssvRf369RWfS/aUnZ0tdu7cKXbu3CkAuOfy5V0Xr7zyiqhTp474+uuvxZ49e8TQoUNF48aNFUvN12rfbreLO+64QyQmJopdu3aJ8+fPux82m0319sui9BqLitr/8ssvhclkEu+99544fPiwePvtt4XBYBAbN27UpP0+ffqINm3aiJ9//lkcO3ZMLF68WAQHB4t//etfirT/6KOPiqioKLFu3TqvP9+8vDz3NY888ohITk4WP/30k/j9999F9+7dRffu3TVpPysrS3Tr1k2kpKSII0eOeF2jxChqZX7+q0HhNRaV6cObb74pIiMjxYoVK8Thw4fFc889J4KDg8WRI0dUb99ut4umTZuKXr16ia1bt4ojR46IuXPnCkmSyl2L4IvJkyeL9evXi+PHj4s9e/aIyZMnC0mSxA8//CCEUPf9V1H7ar//KoPBQgHffPONaNu2rbBYLKJly5bivffe06xtq9UqnnjiCZGcnCyCg4PF9ddfL6ZOnarYh2hZfv75ZwGg1GP06NFCiKItp9OmTRNxcXHCYrGIAQMGiIMHD2rS/vHjx8v8HgDx888/q95+WZQOFpVp/9///rdo2rSpCA4OFu3btxdfffWVZu2fP39ejBkzRiQkJIjg4GDRokUL8frrr7sXs1ZXeX++ixcvdl+Tn58vHnvsMVG3bl0RGhoq7rzzTnH+/HlN2i/v3w8Acfz4cdXbL+81SgaLyvZh9uzZIjExUYSGhoru3bsrFm4r0/6hQ4fE8OHDRWxsrAgNDRXt2rUrtf20qh588EHRqFEjYTabRUxMjBgwYIA7VAih7vuvovbVfv9VBo9NJyIiIsUE9BoLIiIi0hcGCyIiIlIMgwUREREphsGCiIiIFMNgQURERIphsCAiIiLFMFgQERGRYhgsiIiISDEMFkRERKQYBgsiIiJSDIMFERERKYbBgoiIiBTz/01jkXNiR1BMAAAAAElFTkSuQmCC"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_condition_number():\n",
    "    x = np.linspace(4, 32, 15)\n",
    "    y1 = []\n",
    "    y2 = []\n",
    "    for i in range(15):\n",
    "        y1.append(np.linalg.cond(vandermonde_matrix(int(x[i]))))\n",
    "        y2.append(np.linalg.cond(fourier_basis_matrix(int(x[i]))))\n",
    "    print(y1)\n",
    "    print(y2)\n",
    "    plt.plot(x, y1, label='Vandermonde')\n",
    "    plt.plot(x, y2, label='Fourier')\n",
    "    plt.xticks(x)\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "plot_condition_number()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-05-27T09:21:33.202595300Z",
     "start_time": "2023-05-27T09:21:33.026439100Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### explain the reasons for the trends in these two plots\n",
    "As N increases, the condition number of Vandermonde matrix increases, and the condition number of Fourier basis matrix decreases. The reason is that the Vandermonde matrix is ill-conditioned, and the Fourier basis matrix is well-conditioned."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### 3. Construct $A_V = V^T V$ and $A_F = F^T F$ for $N$ = 4, 6, . . . 32, and judge whether they are positive definite or not."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N = 4:(True, True)\n",
      "N = 6:(True, True)\n",
      "N = 8:(True, True)\n",
      "N = 10:(True, True)\n",
      "N = 12:(True, True)\n",
      "N = 14:(False, True)\n",
      "N = 16:(False, True)\n",
      "N = 18:(False, True)\n",
      "N = 20:(False, True)\n",
      "N = 22:(False, False)\n",
      "N = 24:(False, True)\n",
      "N = 26:(False, False)\n",
      "N = 28:(False, False)\n",
      "N = 30:(False, False)\n",
      "N = 32:(False, False)\n"
     ]
    }
   ],
   "source": [
    "def is_positive_definite(n):\n",
    "    vm = vandermonde_matrix(n)\n",
    "    fm = fourier_basis_matrix(n)\n",
    "    vm_t = np.transpose(vm)\n",
    "    fm_t = np.transpose(fm)\n",
    "    vm_a = np.dot(vm_t, vm)\n",
    "    fm_a = np.dot(fm_t, fm)\n",
    "    vm_r = np.all(np.linalg.eigvals(vm_a) > 0)\n",
    "    fm_r = np.all(np.linalg.eigvals(fm_a) > 0)\n",
    "    return vm_r, fm_r\n",
    "\n",
    "for i in range(4, 33, 2):\n",
    "    print(\"N = \" + str(i) + \":\" + str(is_positive_definite(i)))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-05-27T09:36:27.000187100Z",
     "start_time": "2023-05-27T09:36:26.969577100Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### A table of values that includes the following columns: $N, isposdef(A_V ), isposdef(A_F ), cond(V ), cond(F)$\n",
    "\n",
    "| N  | isposdef(A_V) | isposdef(A_F) | cond(V)  | cond(F)  |\n",
    "|----|---------------|---------------|----------|----------|\n",
    "| 4  | True          | True          | 163.4    | 5.027    |\n",
    "| 6  | True          | True          | 8.215e3  | 31.64    |\n",
    "| 8  | True          | True          | 4.5e5    | 236.6    |\n",
    "| 10 | True          | True          | 2.568e7  | 1.918e3  |\n",
    "| 12 | True          | True          | 1.5e9    | 1.631e4  |\n",
    "| 14 | False         | True          | 8.906e10 | 1.43e5   |\n",
    "| 16 | False         | True          | 5.345e12 | 1.283e6  |\n",
    "| 18 | False         | True          | 3.232e14 | 1.171e7  |\n",
    "| 20 | False         | True          | 1.95e16  | 1.083e8  |\n",
    "| 22 | False         | False         | 1.621e18 | 1.012e9  |\n",
    "| 24 | False         | True          | 1.786e19 | 9.536e9  |\n",
    "| 26 | False         | False         | 6.359e18 | 9.055e10 |\n",
    "| 28 | False         | False         | 8.24e19  | 8.651e11 |\n",
    "| 30 | False         | False         | 7.557e18 | 8.309e12 |\n",
    "| 32 | False         | False         | 4.589e18 | 8.049e13 |\n",
    "\n",
    "##### Explain the reasons for the trends in the table\n",
    "\n",
    "The largest value of $N$ where $A_V$ is positive definite: 12, the condition number of that $V$ : 1.5e9\n",
    "The largest value of $N$ where $A_F$ is positive definite: 24, the condition number of that $F$ : 9.536e9\n",
    "These condition numbers connected to the positive definiteness of the matrices. The larger the condition number, the more ill-conditioned the matrix is, and the more likely it is to be not positive definite. This is because the condition number is the ratio of the largest to the smallest singular value of the matrix. The larger the condition number, the more the matrix amplifies the error in the solution.\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### 4. For N = 8, transform the linear systems above into symmetric positive definite systems according to Question 3.3, and use Cholesky factorization to solve them. Report the residual L2 norm for each solution."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2.0202326468819256e-11, 2.0806242997996403e-14)\n"
     ]
    }
   ],
   "source": [
    "# transform the linear systems into symmetric positive definite systems and use Cholesky factorization to solve them\n",
    "def cholesky_factorization(n):\n",
    "    # f(x) = 1 / (1 + x^2)\n",
    "    x = np.linspace(0, 1, n + 1)[:-1]\n",
    "    y = 1 / (1 + x**2)\n",
    "\n",
    "    vm = vandermonde_matrix(n)\n",
    "    fm = fourier_basis_matrix(n)\n",
    "    vm_t = np.transpose(vm)\n",
    "    fm_t = np.transpose(fm)\n",
    "    vm_a = np.dot(vm_t, vm)\n",
    "    fm_a = np.dot(fm_t, fm)\n",
    "    vm_b = np.dot(vm_t, y)\n",
    "    fm_b = np.dot(fm_t, y)\n",
    "    vm_c = sp.linalg.cho_factor(vm_a)\n",
    "    fm_c = sp.linalg.cho_factor(fm_a)\n",
    "    vm_x = sp.linalg.cho_solve(vm_c, vm_b)\n",
    "    fm_x = sp.linalg.cho_solve(fm_c, fm_b)\n",
    "    vm_y = np.dot(vm, vm_x)\n",
    "    fm_y = np.dot(fm, fm_x)\n",
    "    vm_r = np.linalg.norm(vm_y - y)\n",
    "    fm_r = np.linalg.norm(fm_y - y)\n",
    "    return vm_r, fm_r\n",
    "print(cholesky_factorization(8))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-05-27T09:53:55.640220900Z",
     "start_time": "2023-05-27T09:53:55.625011800Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Compare the residuals to Question 3.1: how does Cholesky compare to LU?\n",
    "The residual of Cholesky is bigger than LU. The reason is that the condition number of the Vandermonde matrix is very large, and the Cholesky factorization is not stable for ill-conditioned matrices."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Least Squares Problems and QR"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### 1. Solve the least square system with QR decomposition when M = 16, N = 4, 8."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.58215047 -1.06220587 -0.01392258  1.00103892]\n",
      "[ 1.21467707 -0.11785824  0.35830956  0.48758784]\n",
      "(0.00408323826729538, 0.30230113370126455)\n",
      "[-3.12130078e-02  4.78868165e-01 -1.49113519e+00  1.69087565e+00\n",
      " -1.65663052e-01 -9.80828040e-01 -8.07368955e-04  1.00000139e+00]\n",
      "[ 1.33778075 -0.1950913  -0.31748875  0.02996036  0.38439885  0.78902919\n",
      " -0.10683444 -0.06818281]\n",
      "(1.341856782350059e-05, 0.008863804298810679)\n"
     ]
    }
   ],
   "source": [
    "# construct a m*n vandermonde matrix\n",
    "def vandermonde_matrix_mn(m,n):\n",
    "    x = np.linspace(0, 1, m + 1)[:-1]\n",
    "    return np.vander(x, n)\n",
    "def fourier_basis_matrix_mn(m, n):\n",
    "    x = np.linspace(0, 1, m + 1)[:-1]\n",
    "    y = np.zeros((m, n))\n",
    "    for i in range(m):\n",
    "        for j in range(n):\n",
    "            if j + 1 <= n/2:\n",
    "                y[i][j] = np.sin((j + 1) * np.pi * x[i])\n",
    "            else:\n",
    "                y[i][j] = np.cos((j + 1 - n/2) * np.pi * x[i])\n",
    "    return y\n",
    "# solve the least square system with QR decomposition by minimizing the 2-norm of the residual\n",
    "def least_squares_qr(m, n):\n",
    "    x = np.linspace(0, 1, m + 1)[:-1]\n",
    "    y = 1 / (1 + x**2)\n",
    "    vm = vandermonde_matrix_mn(m, n)\n",
    "    fm = fourier_basis_matrix_mn(m, n)\n",
    "    vm_q, vm_r = np.linalg.qr(vm)\n",
    "    fm_q, fm_r = np.linalg.qr(fm)\n",
    "    vm_qt = np.transpose(vm_q)\n",
    "    fm_qt = np.transpose(fm_q)\n",
    "    vm_b = np.dot(vm_qt, y)\n",
    "    fm_b = np.dot(fm_qt, y)\n",
    "    vm_x = sp.linalg.solve_triangular(vm_r, vm_b)\n",
    "    fm_x = sp.linalg.solve_triangular(fm_r, fm_b)\n",
    "    print(vm_x)\n",
    "    print(fm_x)\n",
    "    vm_y = np.dot(vm, vm_x)\n",
    "    fm_y = np.dot(fm, fm_x)\n",
    "    vm_r = np.linalg.norm(vm_y - y)\n",
    "    fm_r = np.linalg.norm(fm_y - y)\n",
    "    return vm_r, fm_r\n",
    "print(least_squares_qr(16, 4))\n",
    "print(least_squares_qr(16, 8))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-05-27T12:06:46.143005900Z",
     "start_time": "2023-05-27T12:06:46.034646Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### 2. Plot the $g_V$ , $g_F$ when M = 16, N = 4, 8, compare them with the analytical function f(x) and the interpolation function obtained in Question 3.1."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGzCAYAAAD9pBdvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABwpklEQVR4nO3deVhUZfsH8O8wMAMjmwYMYCjgjpoLiispSZKa2m654r6UGz9NzV5J63VJQ8slytxSy90yNU0oCs1eTEQzhTQR0AQ02QRlcOb5/UFMjgM4gwNH4Pu5rrl0nnnOOfecmWHuebYjE0IIEBEREUnESuoAiIiIqHZjMkJERESSYjJCREREkmIyQkRERJJiMkJERESSYjJCREREkmIyQkRERJJiMkJERESSYjJCREREkmIyQjXGO++8A5lMZlbdGzduVHJU1VNMTAxkMhliYmKkDqVcoaGh8Pb2ljqMasvb2xuhoaEm1e3Zsyd69uxZqfFQ7VWtkpGNGzdCJpPh119/lToUA+fOncM777yDy5cvm7zN0aNH0adPH9SvXx+2trZo0KAB+vfvjy+++KLyAq2FFi5ciK+++qrS9r9//34888wzeOyxx2Bra4umTZti5syZuHnzplHd0NBQyGQy/U2pVKJp06aYN28e7ty5U2kxVqYHfSYvX74MmUyGZcuWlfr4smXLIJPJzPrsUOWpyN+yylaSGMtkMmzZsqXUOt26dYNMJkOrVq2qNLZr165h9uzZCAoKgoODwwMTeI1Gg4ULF6J58+awtbWFWq1Gv379cOXKlXKPU/I5kslk2L17t9Hjj8KPq61bt0Imk8He3r5C21erZORRde7cOcyfP9/kD/DOnTvx5JNPIiMjA1OnTsXKlSsxdOhQZGVlYe3atZUbbA329ttv4/bt2wZllZmMzJgxA/3790d6ejpmzZqFVatWITg4GCtXrkTbtm1x4cIFo22USiU2b96MzZs3IyIiAt7e3nj33XcxevToSomRqDxJSUkGf3PK+1v23Xff4bvvvqvC6AzZ2tqW+mPt8uXL+Pnnn2Fra1vlMSUlJWHJkiW4evUqWrduXW7doqIi9OvXD//973/xzDPPYM2aNXjzzTdRp04d5OTkmHzMBQsW4FG7pNytW7f0z6WirC0YD5nonXfegZ+fH3755RcoFAqDxzIzMyWK6sF0Oh00Go0kH3pTWFtbw9q6at7SX375JT744AMMGjQIW7duhVwu1z8WGhqKoKAgvPzyy/j1118NYrK2tsbQoUP19ydNmoSuXbviyy+/REREBNRqdZXETwQUJ8emuv9vVVXr27cv9u3bhxs3bsDFxUVf/sUXX0CtVqNJkybIysqq0pj8/f3x999/o169eti1axdefvnlMusuX74cP/74I44ePYqAgIAKHa9t27ZISEjA3r178cILL1Q0bIt777334ODggKCgoAr/+KuRLSNXr17FqFGjoFaroVQq0bJlS6xfv96gjkajwbx58+Dv7w8nJyfUqVMHgYGB+OGHH4z2t23bNvj7+8PBwQGOjo5o3bo1PvzwQwDFzdQlb8CgoCB9U1p5TXV//vknOnbsWOqH283NzeB+dnY2QkND4eTkBGdnZ4wYMQIJCQmQyWTYuHGjvl5Z/bml9akvW7YMXbt2xWOPPQY7Ozv4+/tj165dRtvKZDK88cYb2Lp1K1q2bAmlUolDhw4BMO0cA8DKlSvRsmVLqFQq1K1bFx06dCi3K0oIARcXF4SFhenLdDodnJ2dIZfLkZ2drS9fsmQJrK2tcevWLQDGY0ZkMhny8/OxadMm/etyf/94yfl1dnaGk5MTRo4ciYKCgjLjKzF//nzUrVsXn376qUEiAgABAQGYNWsWTp8+jT179pS7H5lMhu7du0MIgUuXLj3wuF9//TX69esHT09PKJVKNGrUCO+++y60Wq1BvZ49e6JVq1Y4d+4cgoKCoFKpUL9+fbz//vtG+7xy5Qqee+451KlTB25ubpg+fToKCwsfGEtV++qrr9CqVSvY2tqiVatW2Lt3b6n1dDodVqxYgZYtW+qbwsePH1/qF9W3336LHj166D/bHTt2NHp/7ty5E/7+/rCzs4OLiwuGDh2Kq1evGtQJDQ2Fvb09UlNT8eyzz8Le3h7169fH6tWrAQC//fYbnnrqKdSpUwcNGzY0OkZJd9dPP/2E8ePH47HHHoOjoyOGDx9eatxr1qzRfyY9PT3x+uuvG3w2AODChQt48cUX4e7uDltbWzz++ON49dVXDX6F3ztm5EF/y0r7G5OZmYnRo0dDrVbD1tYWbdq0waZNmwzq3NtV9+mnn6JRo0ZQKpXo2LEjTpw4YfTcyjJw4EAolUrs3LnToPyLL77AK6+8YvQ5rAoODg6oV6/eA+vpdDp8+OGHeP755xEQEIC7d++a9Hfmfq+++iqaNm36SLWOXLhwAcuXL0dERMRD/RiscS0jGRkZ6Ny5s/6L1NXVFd9++y1Gjx6N3NxcTJs2DQCQm5uLzz77DK+99hrGjh2LvLw8rFu3DiEhIYiLi0Pbtm0BAEeOHMFrr72GXr16YcmSJQCA8+fP49ixY5g6dSqefPJJTJkyBR999BHeeusttGjRAgD0/5amYcOGiI6OxpUrV/D444+XWU8IgYEDB+Lo0aOYMGECWrRogb1792LEiBEPdY4+/PBDDBgwAEOGDIFGo8G2bdvw8ssvY//+/ejXr59B3e+//x47duzAG2+8ARcXF3h7e5t8jteuXYspU6bgpZdewtSpU3Hnzh2cOXMG//vf/zB48OBSY5PJZOjWrRt++uknfdmZM2eQk5MDKysrHDt2TB9jbGws2rVrV2Yf5ebNmzFmzBgEBARg3LhxAIBGjRoZ1HnllVfg4+ODRYsWIT4+Hp999hnc3Nz0r3VpLly4gKSkJISGhsLR0bHUOsOHD0d4eDi++eYbvPLKK2XuC4C+Sbxu3brl1gOKvzDs7e0RFhYGe3t7fP/995g3bx5yc3OxdOlSg7pZWVl45pln8MILL+CVV17Brl27MGvWLLRu3Rp9+vQBANy+fRu9evVCamoqpkyZAk9PT2zevBnff//9A2OpSt999x1efPFF+Pn5YdGiRfj7778xcuTIUj8/48ePx8aNGzFy5EhMmTIFycnJWLVqFU6dOoVjx47BxsYGQPG5HDVqFFq2bIk5c+bA2dkZp06dwqFDh/Tvz5L9dOzYEYsWLUJGRgY+/PBDHDt2DKdOnYKzs7P+uFqtFn369MGTTz6J999/H1u3bsUbb7yBOnXqYO7cuRgyZAheeOEFREZGYvjw4ejSpQt8fHwMYn/jjTfg7OyMd955B0lJSfj444+RkpKiHzcBFCfd8+fPR3BwMCZOnKivd+LECf3z02g0CAkJQWFhISZPngx3d3dcvXoV+/fvR3Z2NpycnIzOm7l/y27fvo2ePXvi4sWLeOONN+Dj44OdO3ciNDQU2dnZmDp1qkH9L774Anl5eRg/fjxkMhnef/99vPDCC7h06ZL+NSmPSqXCwIED8eWXX2LixIkAgNOnT+P333/HZ599hjNnzjxwHwBQUFBgUiIgl8tN+kya4ty5c/jrr7/wxBNPYNy4cdi0aRM0Go3+h21QUJBJ+5HL5Xj77bcxfPjwCrWOFBUVmdwlVK9ePVhZPbi9Ytq0aQgKCkLfvn2xY8cOs+IxIKqRDRs2CADixIkTZdYZPXq08PDwEDdu3DAof/XVV4WTk5MoKCgQQghx9+5dUVhYaFAnKytLqNVqMWrUKH3Z1KlThaOjo7h7926Zx9y5c6cAIH744QeTnse6desEAKFQKERQUJD4z3/+I2JjY4VWqzWo99VXXwkA4v3339eX3b17VwQGBgoAYsOGDfryHj16iB49ehgda8SIEaJhw4YGZSXnoIRGoxGtWrUSTz31lEE5AGFlZSV+//13g3JTz/HAgQNFy5Ytyz0XpVm6dKmQy+UiNzdXCCHERx99JBo2bCgCAgLErFmzhBBCaLVa4ezsLKZPn67fLjw8XNz/lq5Tp44YMWKE0TFK6t77WgshxPPPPy8ee+yxcuMreV2WL19ebj1HR0fRvn17/f0RI0aIOnXqiOvXr4vr16+LixcvimXLlgmZTCZatWoldDpdufsTwvi1E0KI8ePHC5VKJe7cuaMv69GjhwAgPv/8c31ZYWGhcHd3Fy+++KK+bMWKFQKA2LFjh74sPz9fNG7c2KT39IM+k8nJyQKAWLp0aamPL126VAAQycnJ5R6nbdu2wsPDQ2RnZ+vLvvvuOwHA4P0dGxsrAIitW7cabH/o0CGD8uzsbOHg4CA6deokbt++bVC35HXQaDTCzc1NtGrVyqDO/v37BQAxb948fdmIESMEALFw4UJ9WVZWlrCzsxMymUxs27ZNX56YmCgAiPDwcH1ZyXn09/cXGo1GX/7+++8LAOLrr78WQgiRmZkpFAqF6N27t8Hfi1WrVgkAYv369UIIIU6dOiUAiJ07d5Z7Xhs2bGjw+Sjvb9n9f2NK3jtbtmzRl2k0GtGlSxdhb2+v//yWvAcee+wxcfPmTX3dr7/+WgAQ33zzTbkx/vDDD/rnsn//fiGTyURqaqoQQoiZM2cKX19ffXym/L0p+ew/6Hb/380HKe/c7dmzR38OmjRpIjZs2CA2bNggmjRpIhQKhTh9+nS5+773c3T37l3RpEkT0aZNG/17teQ5Xb9+vdz9lJxLU24P+kwKUfxZsLa21n9HlPyNq4ga1U0jhMDu3bvRv39/CCFw48YN/S0kJAQ5OTmIj48HUJxhlnST6HQ63Lx5E3fv3kWHDh30dQDA2dkZ+fn5OHLkiMXiHDVqFA4dOoSePXvi6NGjePfddxEYGIgmTZrg559/1tc7ePAgrK2t9b8CSuKePHnyQx3fzs5O//+srCzk5OQgMDDQ4HmX6NGjB/z8/PT3zTnHzs7OuHLlillNsQAQGBgIrVarPxexsbEIDAxEYGAgYmNjAQBnz55FdnY2AgMDzX7+95owYYLRsf/++2/k5uaWuU1eXh6A4iba8jg4OOjrlsjPz4erqytcXV3RuHFjzJgxA926dcPXX39t0rTke1+7vLw83LhxA4GBgSgoKEBiYqJBXXt7e4PxKQqFAgEBAQbdQQcPHoSHhwdeeuklfZlKpdK3JD0Krl27hoSEBIwYMcLgF/3TTz9t8N4EirtUnJyc8PTTTxu8N/39/WFvb6/vhj1y5Ajy8vIwe/ZsozFQJa/Dr7/+iszMTEyaNMmgTr9+/dC8eXMcOHDAKNYxY8bo/+/s7IxmzZqhTp06Bq1jzZo1g7Ozc6ndcuPGjTNoJZg4cSKsra1x8OBBAEBUVBQ0Gg2mTZtm8Kt17NixcHR01MdUcp4OHz5coe4AUxw8eBDu7u547bXX9GU2NjaYMmUKbt26hR9//NGg/qBBgwxaGko+u6Z0T5bo3bs36tWrh23btkEIgW3bthkc3xTDhw/HkSNHHnjbunWrWfstT0lXcl5eHqKjoxEaGorQ0FBERUVBCFFq92lZSlpHTp8+bfb4jDZt2pj03I8cOQJ3d/dy96XRaDB9+nRMmDDB6HNYETWqm+b69evIzs7Gp59+ik8//bTUOvcOEN20aRM++OADJCYmoqioSF9+b9PppEmTsGPHDv003N69e+OVV17BM88881CxhoSEICQkBAUFBTh58iS2b9+OyMhIPPvss0hMTISbmxtSUlLg4eFh1A3RrFmzhzr2/v378d577yEhIcFgbEBpX4b3NyObc45nzZqFqKgoBAQEoHHjxujduzcGDx6Mbt26lRtf+/btoVKpEBsbi5CQEMTGxmL+/Plwd3fHypUrcefOHX1S0r17d7Oe+/0aNGhgcL/kj2VWVlaZXTAlScj9icb98vLyjMbr2Nra4ptvvgFQPFbj/fffR2ZmpkGSUZ7ff/8db7/9Nr7//nujhOn+5tfHH3/c6DWtW7euQXN2SkoKGjdubFTvYd9j5iovEUtJSQEANGnSxOixZs2aGSTRFy5cQE5OjtHYqxIl780///wTAMqdClpy3NLORfPmzXH06FGDMltbW7i6uhqUOTk5lfo6ODk5lToW5P7naG9vDw8PD31XXlkxKRQK+Pr66h/38fFBWFgYIiIisHXrVgQGBmLAgAEYOnRoqV00FZGSkoImTZoYNeWXdOuUxFKivM+aqWxsbPDyyy/jiy++QEBAANLS0srs8i2Lr68vfH19zdrmYZV8vrt16wYvLy99eYMGDdC9e3eDH6GmGDJkCN59910sWLAAzz33nMnb1a1bF8HBwWYdqyzLly/HjRs3MH/+fIvsr0YlIzqdDgAwdOjQMsdVPPHEEwCALVu2IDQ0FM899xxmzpwJNzc3yOVyLFq0SP+HCigeUJqQkIDDhw/j22+/xbfffosNGzZg+PDhRgO1KkKlUul/9bu4uGD+/Pn49ttvzR4XIpPJSh3QdP/AxtjYWAwYMABPPvkk1qxZAw8PD9jY2GDDhg2lDiy9/0vSnHPcokULJCUlYf/+/Th06BB2796NNWvWYN68eeW+gW1sbNCpUyf89NNPuHjxItLT0xEYGAi1Wo2ioiL873//Q2xsLJo3b270x99cZQ16K+1clij5FVBeH3VKSgpyc3ON/ujJ5XKDPwYhISFo3rw5xo8fj3379pUba3Z2Nnr06AFHR0csWLAAjRo1gq2tLeLj4zFr1iz9a/Mwz83SSloU7p9yXaLkV7ulZmjpdDq4ubmV+av2Yd8v5SnrfEv1OnzwwQcIDQ3F119/je+++w5TpkzBokWL8Msvv5Q7Vq2yWOo8DB48GJGRkXjnnXfQpk0bs3+V37p1S99SUR65XG6x94unpycAlDpbzs3NDadOnTJrfyWtIyWvr6k0Gk2payCVxtXVtczXLCcnB++99x4mTZqE3Nxc/Q+jW7duQQiBy5cvQ6VSlfmjoDQ1KhlxdXWFg4MDtFrtA7O/Xbt2wdfXF3v27DH41RIeHm5UV6FQoH///ujfvz90Oh0mTZqETz75BP/5z39K/VVZUR06dABQ3CwN/DvQ9datWwatI0lJSUbb1q1bt9Tmzvt/nezevRu2trY4fPiwwbS+DRs2mBSjOecYAOrUqYNBgwZh0KBB0Gg0eOGFF/Df//4Xc+bMKfcLKDAwEEuWLEFUVBRcXFzQvHlzyGQytGzZErGxsYiNjcWzzz77wONb6rW5V5MmTdCsWTN89dVX+PDDD0vtrvn8888BoNypfgDg4eGB6dOnY/78+fjll1/QuXPnMuvGxMTg77//xp49e/Dkk0/qy5OTkyv4TIrfY2fPnoUQwuBclfYeqwhXV1eoVKoy95eUlASVSmUwVbO0GAGUum7L/ftt1KgRoqKi0K1bt3Jbm0oGMp89exaNGzcu97hJSUl46qmnjI5b8rglXbhwwWAw461bt3Dt2jX07dvXKKZ7E12NRoPk5GSjz2Tr1q3RunVrvP322/j555/RrVs3REZG4r333iv1+OZ8Xho2bIgzZ85Ap9MZtI6UdBdWxvkBiltDGzRogJiYmHIHmpdl2bJlJv2ab9iwocUWf2vdujVsbGyMZmEBwF9//VWhpGfo0KF47733MH/+fAwYMMCkbX7++WeTB8smJyeXubpxVlYWbt26hffff7/ULiYfHx8MHDjQrG6kGjVmRC6X48UXX8Tu3btx9uxZo8evX79uUBcwzMr/97//4fjx4wbb/P333wb3rays9L/8S7o4ShZ6uX9qXVmio6NLLS/pFy5pgu3bty/u3r2Ljz/+WF9Hq9Vi5cqVRts2atQIiYmJBs/x9OnTOHbsmEE9uVwOmUxm0GJy+fJlk9805pzj+8+dQqGAn58fhBAG3WKlCQwMRGFhIVasWIHu3bvr/0gGBgZi8+bN+Ouvv0waL1KnTh2TXxdzhIeHIysrCxMmTDBqfTp58iSWLFmCdu3a6WetlGfy5MlQqVRYvHhxufVKe89qNBqsWbOmAs+gWN++ffHXX38ZTO0uKCgoswvOXHK5HL1798Y333yD1NRUg8dSU1PxzTffoHfv3uVOy/Tw8EDbtm2xadMmg66oI0eO4Ny5cwZ1X3nlFWi1Wrz77rtG+7l7967+vdC7d284ODhg0aJFRqvflpzfDh06wM3NDZGRkQbdmd9++y3Onz9vNPPMEj799FODz8bHH3+Mu3fv6t9HwcHBUCgU+OijjwzeB+vWrUNOTo4+ptzcXNy9e9dg361bt4aVlVW507bN+VvWt29fpKenY/v27fqyu3fvYuXKlbC3t0ePHj0e/IQrQCaT4aOPPkJ4eDiGDRtm9vZSjBlxcHBA37598fPPPxuM7Tp//jx+/vlnPP3002bvs6R1JCEh4YGtqiUsNWbEzc0Ne/fuNboFBQXB1tYWe/fuxZw5c8x6PtWyZWT9+vX69S7uNXXqVCxevBg//PADOnXqhLFjx8LPzw83b95EfHw8oqKi9E1Uzz77LPbs2YPnn38e/fr1Q3JyMiIjI+Hn52fQhDdmzBjcvHkTTz31FB5//HGkpKToV9gs6Rtt27Yt5HI5lixZgpycHCiVSjz11FNlNlENHDgQPj4+6N+/Pxo1aoT8/HxERUXhm2++QceOHdG/f38AQP/+/dGtWzfMnj0bly9fhp+fH/bs2VPq1KxRo0YhIiICISEhGD16NDIzMxEZGYmWLVsajC3o168fIiIi8Mwzz2Dw4MHIzMzE6tWr0bhxY5Onxpl6jnv37g13d3d069YNarUa58+fx6pVq9CvX78HDv7s0qULrK2tkZSUZDCY8sknn9QnZ6YkI/7+/oiKikJERAQ8PT3h4+ODTp06mfQ8y/Paa6/h119/RUREBM6dO4chQ4agbt26iI+Px/r16+Hq6opdu3aZNO/+sccew8iRI7FmzRqcP3++zKmUXbt2Rd26dTFixAhMmTIFMpkMmzdvfqjm/rFjx2LVqlUYPnw4Tp48CQ8PD2zevBkqlcqs/ZT3mVy4cCE6d+6M9u3bY9y4cfD29sbly5fx6aefQiaTYeHChQ/c/6JFi9CvXz90794do0aNws2bN/Vr2Nz7ee3RowfGjx+PRYsWISEhAb1794aNjQ0uXLiAnTt34sMPP8RLL70ER0dHLF++HGPGjEHHjh0xePBg1K1bF6dPn0ZBQQE2bdoEGxsbLFmyBCNHjkSPHj3w2muv6af2ent7Y/r06WadI1NoNBr06tULr7zyCpKSkrBmzRp0795d/8vX1dUVc+bMwfz58/HMM89gwIAB+nodO3bUD1j+/vvv8cYbb+Dll19G06ZNcffuXWzevFn/Y6Is5vwtGzduHD755BOEhobi5MmT8Pb2xq5du3Ds2DGsWLHigZ/xhzFw4EAMHDiwQttaesxISSvT77//DqB4SYGS8URvv/22vt7ChQsRHR2Np556ClOmTAEAfPTRR6hXrx7eeuutCh27ZOxIQkKCSfUtNWZEpVKVOlblq6++QlxcnFnjWPQqNAdHIiXT38q6paWlCSGEyMjIEK+//rrw8vISNjY2wt3dXfTq1Ut8+umn+n3pdDqxcOFC0bBhQ6FUKkW7du3E/v37jabC7tq1S/Tu3Vu4ubkJhUIhGjRoIMaPHy+uXbtmENvatWuFr6+vkMvlD5wS+eWXX4pXX31VNGrUSNjZ2QlbW1vh5+cn5s6dq58OV+Lvv/8Ww4YNE46OjsLJyUkMGzZMP23v3qm9QgixZcsW4evrKxQKhWjbtq04fPhwqVN7161bJ5o0aSKUSqVo3ry52LBhQ6nTYgGI119/vdTnYMo5/uSTT8STTz4pHnvsMaFUKkWjRo3EzJkzRU5OTpnn5l4dO3YUAMT//vc/fdmVK1cEAOHl5WVUv7TnkJiYKJ588klhZ2cnAOinMZY1Fa7kPWbKtDYhhNi3b58IDg4Wzs7O+vdhy5YtS32O5U17+/PPP4VcLi91GvK9jh07Jjp37izs7OyEp6enePPNN8Xhw4eN3nNlTXMs7f2QkpIiBgwYIFQqlXBxcRFTp07VT4U1dWrvgz6T58+fF4MGDRJubm7C2tpauLm5iVdffVWcP3++3P3fa/fu3aJFixZCqVQKPz8/sWfPnlKfjxBCfPrpp8Lf31/Y2dkJBwcH0bp1a/Hmm2+Kv/76y6Devn37RNeuXYWdnZ1wdHQUAQEB4ssvvzSos337dtGuXTuhVCpFvXr1xJAhQ8SVK1cM6pT12pb1OjRs2FD069fP6Dz++OOPYty4caJu3brC3t5eDBkyRPz9999G269atUo0b95c2NjYCLVaLSZOnCiysrL0j1+6dEmMGjVKNGrUSNja2op69eqJoKAgERUVZRTH/e+5sv6WlbZ8QEZGhhg5cqRwcXERCoVCtG7d2ujvUnnTu3HfFOfS3Du1tzymTu21tPLe//c7efKkCA4OFnXq1BEODg5i4MCB4o8//njgMco7h/d+Bh80tbeyPczUXpkQj8gybmSyy5cvw8fHBxs2bDD5iptUNcaMGYN169Zh7dq1BtM8icpTsrjaiRMn9GPHiGqTatlNQ/So+uSTT5CRkYGJEyfC09NTP/CQiIjKxmSEyILkcrl+HREiIjJNjZpNQ0RERNUPx4wQERGRpNgyQkRERJJiMkJERESSqhYDWHU6Hf766y84ODhUyvLeREREZHlCCOTl5cHT09Pooor3qhbJyF9//WVwpUMiIiKqPtLS0sq9QGO1SEZKlhVOS0sr87LuRERE9GjJzc2Fl5fXAy8PUC2SkZKuGUdHRyYjRERE1cyDhlhwACsRERFJiskIERERSYrJCBEREUmqWowZISKimkkIgbt370Kr1UodClWAXC6HtbX1Qy+7wWSEiIgkodFocO3aNRQUFEgdCj0ElUoFDw8PKBSKCu+DyQgREVU5nU6H5ORkyOVyeHp6QqFQcFHLakYIAY1Gg+vXryM5ORlNmjQpd2Gz8jAZISKiKqfRaKDT6eDl5QWVSiV1OFRBdnZ2sLGxQUpKCjQaDWxtbSu0Hw5gJSIiyVT0lzQ9OizxGrJlpJbQ6rSIz4zH9YLrcFW5or1be8it5FKHRURExGSkNohKicLiuMXIKMjQl6lVaswOmI3ghsESRkZERMRumhovKiUKYTFhBokIAGQWZCIsJgxRKVESRUZE9PC0OoHjf/6NrxOu4viff0OrE5V+TCEExo0bh3r16kEmkyEhIQF///033NzccPnyZZP2odFo4O3tjV9//bVyg60m2DJSg2l1WiyOWwwB4w+ngIAMMiyJW4IgryB22RBRtXPo7DXM/+YcruXc0Zd5ONkivL8fnmnlUXnHPXQIGzduRExMDHx9feHi4oI333wTAwcOhLe3t0n7UCgUmDFjBmbNmoXo6OhKi7W6YMtIDRafGW/UInIvAYH0gnTEZ8ZXWgxanRYn0k/g4KWDOJF+AlodFzYiood36Ow1TNwSb5CIAEB6zh1M3BKPQ2evVdqx//zzT3h4eKBr165wd3eHRqPBunXrMHr0aLP2M2TIEBw9ehS///57JUVafTAZqcGuF1y3aD1zRaVEIWR3CEYdHoVZsbMw6vAohOwOYdcQET0UrU5g/jfnSmnzhb5s/jfnKqXLJjQ0FJMnT0ZqaipkMhm8vb1x8OBBKJVKdO7cWV9vwYIF8PT0xN9//60v69evH4KCgqDT6QAAdevWRbdu3bBt2zaLx1ndMBmpwVxVrhatZw6OVSGiyhKXfNOoReReAsC1nDuIS75p8WN/+OGHWLBgAR5//HFcu3YNJ06cQGxsLPz9/Q3qzZ07F97e3hgzZgwAYPXq1fj555+xadMmg6mwAQEBiI2NtXic1Q2TkRqsvVt7qFVqyFD6qoYyyOCuckd7t/YWPe6DxqoAwJK4JeyyIaIKycwrOxGpSD1zODk5wcHBAXK5HO7u7nB1dUVKSgo8PT0N6snlcmzZsgXR0dGYPXs2Zs6cidWrV6NBgwYG9Tw9PZGSkmLxOKsbJiM1mNxKjtkBswHAKCEpuT8rYJbFB68+CmNViKjmcnMwbZVPU+s9rNu3b5e68qivry+WLVuGJUuWYMCAARg8eLBRHTs7O16bB0xGarzghsGI6BkBN5WbQblapUZEz4hKWWdE6rEqRFSzBfjUg4eTbRltvoAMxbNqAnzqVUk8Li4uyMrKKvWxn376CXK5HJcvX8bdu3eNHr958yZcXS3fVV7dMBmpBYIbBuPwi4exPmQ9lgQuwfqQ9Tj04qFKW/BMyrEqRFTzya1kCO/vBwBGCUnJ/fD+fpBbVc2F99q1a4dz584ZlW/fvh179uxBTEwMUlNT8e677xrVOXv2LNq1a1cVYT7SmIzUEnIrOTq6d0Rf377o6N6xUtcVkWqsChHVHs+08sDHQ9vD3cmwe8TdyRYfD21fqeuM3C8kJAS///67QevIlStXMHHiRCxZsgTdu3fHhg0bsHDhQvzyyy8G28bGxqJ3795VFuujiouekcWVjFUJiwmDDDKDgayVOVaFiGqXZ1p54Gk/d8Ql30Rm3h24ORR3zVRVi0iJ1q1bo3379tixYwfGjx8PIQRCQ0MREBCAN954A0BxwjJx4kQMHToUCQkJsLe3x/Hjx5GTk4OXXnqpSuN9FMmEEJW/du5Dys3NhZOTE3JycuDo6Ch1OGSi0q6J465yx6yAWbwmDlEtd+fOHSQnJ8PHx6fCl51/lBw4cAAzZ87E2bNnTb6K7aBBg9CmTRu89dZblRxd5SrvtTT1+5stI1RpghsGI8griFcLJqIar1+/frhw4QKuXr0KLy+vB9bXaDRo3bo1pk+fXgXRPfrYMkJERFWuprWM1GaWaBnhAFYiIiKSFJMRIiIikhSTESIiIpIUkxEiIiKSFJMRIiIikhSTESIiIpIUkxEiIiKSFJMRIiKqvnRaIDkW+G1X8b86rdQRmWXYsGFYuHChyfVnz56NyZMnV2JE0mAyQkRE1dO5fcCKVsCmZ4Hdo4v/XdGquFwirVu3xoQJE0p9bPPmzVAqlbhx4wYA4PTp0zh48CCmTJli8v5nzJiBTZs24dKlS2bFtWnTJnTs2BEqlQoODg7o0aMH9u/fb1AnJiYGMplMf3N1dUXfvn3x22+/mXWsimAyQkRE1c+5fcCO4UDuX4bludeKyyVKSEaPHo1t27bh9u3bRo9t2LABAwYMgIuLCwBg5cqVePnll2Fvb2/y/l1cXBASEoKPP/7Y5G1mzJiB8ePHY9CgQThz5gzi4uLQvXt3DBw4EKtWrTKqn5SUhGvXruHw4cMoLCxEv379oNFoTD5eRTAZISKi6kWnBQ7NAlDa1Uz+KTs0u9K6bPLy8jBkyBDUqVMHHh4eWL58OXr27Ilp06Zh6NChuH37Nnbv3m2wTXJyMmJiYjB69GgAgFarxa5du9C/f399ncTERKhUKnzxxRf6sh07dsDOzg7nzp3Tl/Xv3x/btm0zKdZffvkFH3zwAZYuXYoZM2agcePGaNGiBf773/9i2rRpCAsLQ1pamsE2bm5ucHd3R/v27TFt2jSkpaUhMTHR7PNkDiYjRERUvaT8bNwiYkAAuVeL61WCsLAwHDt2DPv27cORI0cQGxuL+Ph4AMUtFwMHDsT69esNttm4cSMef/xx9O7dGwBw5swZ5OTkoEOHDvo6zZs3x7JlyzBp0iSkpqbiypUrmDBhApYsWQI/Pz99vYCAAFy5cgWXL19+YKxffvkl7O3tMX78eKPH/u///g9FRUVGiVOJnJwcfdKjUCgeeKyHUXuv2qvTFr9Rb2UA9mqgYVeAV5MlInr03cqwbD0z5OXlYdOmTfjiiy/Qq1cvAMXdL56envo6o0ePRp8+ffQXjxNCYNOmTRgxYgSsrIrbAFJSUiCXy+Hm5maw/0mTJuHgwYMYOnQoFAoFOnbsaDRgteRYKSkp8Pb2LjfeP/74A40aNSo1mfD09ISjoyP++OMPg/LHH38cAJCfnw8AGDBgAJo3b/6gU/NQamcycm4f7h6Yjct/uuJWUT3Y29yEd6PrsO63GPAbIHV0RERUHnu1ZeuZ4dKlSygqKkJAQIC+zMnJCc2aNdPff/rpp/H4449jw4YNWLBgAaKjo5GamoqRI0fq69y+fRtKpRIymczoGOvXr0fTpk1hZWWF33//3aiOnZ0dAKCgoMCkmIUorTvrX/cnKrGxsVCpVPjll1+wcOFCREZGmnSch2F2N81PP/2E/v37w9PTEzKZDF999dUDt4mJiUH79u2hVCrRuHFjbNy4sQKhWsi5fTgbsRYbzy3AYYTjmM1kHEY4Np5bgLMRayUdhU1ERCZo2BVw9ARg/EVeTAY41i+uJwErKyuEhoZi06ZN0Ol02LBhA4KCguDr66uv4+LigoKCglIHhp4+fRr5+fnIz8/HtWvXjB6/efMmAMDV1fWBsTRp0gSXLl0q9Th//fUXcnNz0bRpU4NyHx8fNGvWDCNGjMCYMWMwaNCgBx7nYZmdjOTn56NNmzZYvXq1SfWTk5PRr18/BAUFISEhAdOmTcOYMWNw+PBhs4N9aDotEj7aih/l01CoqGvwUKGiLn6UT0PCR1srdZ665vYdHF3+Gb6bswxHl38Gze07lXYsIqIayUoOPLPknzv3JyT/3H9mcaV0vfv6+sLGxgYnTpzQl+Xk5Bh1dYwcORJpaWnYs2cP9u7dqx+4WqJt27YAYDAwFShONEJDQzF37lyEhoZiyJAhRjNzzp49CxsbG7Rs2fKB8b722mu4desWPvnkE6PHli1bBltb23KTjddffx1nz57F3r17H3ish2F2N02fPn3Qp08fk+tHRkbCx8cHH3zwAQCgRYsWOHr0KJYvX46QkBBzD/9QNBd+wq+FLwMKAPc3jclkgBD4tfBl+F34CYpmQRY/fnT4ciSnNkCh8p/sOAtInHQAPg1S0Wv+dIsfr7bT6rSIz4zH9YLrcFW5or1be8g5LoioZvAbALzyefGsmnsHszp6FicildTl7uDggBEjRmDmzJmoV68e3NzcEB4eDisrK4PuFB8fHzz11FMYN24clEolXnjhBYP9uLq6on379jh69Kg+MQGACRMmwMvLC2+//TYKCwvRrl07zJgxw6ABIDY2FoGBgfrumvJ06dIFU6dOxcyZM6HRaPDcc8+hqKgIW7ZswUcffYSNGzfiscceK3N7lUqFsWPHIjw8HM8991yp3UqWUOmzaY4fP47g4GCDspCQEBw/frzMbQoLC5Gbm2tws4STX59GobKecSJSQiZDobIeTn592iLHu1d0+HIkpj+BQoWzQXmhwhmJ6U8gOny5xY9Zm0WlRCFkdwhGHR6FWbGzMOrwKITsDkFUSpTUoRGRpfgNAKadBUbsB15cV/zvtN8qfexfREQEunTpgmeffRbBwcHo1q0bWrRoAVtbW4N6o0ePRlZWFgYPHmz0GACMGTMGW7du1d///PPPcfDgQWzevBnW1taoU6cOtmzZgrVr1+Lbb7/V19u2bRvGjh1rcrwrVqzAmjVr8OWXX6JVq1Zo0aIFli5diu+//x5Dhw594PZvvPEGzp8/j507d5p8THNVejKSnp4OtdpwEJFarUZubm6pi8IAwKJFi+Dk5KS/eXl5WSSW7Fwbi9Yzleb2HVxKbVB8p7QWGQCXUr3YZWMhUSlRCIsJQ0aB4Uj6zIJMhMWEMSEhqkms5IBPIND6peJ/q6D108HBAVu3btWP6Rg3bhySkpLQuHFjg3qvvfYahBBlDmsIDQ3F1atX9T/Ohw8fjlu3bqFJkyb6OgEBAdBoNPoeiW+//RZWVlZ46aWXzIp51KhR+PXXX3H79m0kJyfD3d0da9asgVb777CEnj17QggBZ2dng229vLxQVFSEV155xaxjmuORXGdkzpw5yMnJ0d/uX5ClohTeTR5cyYx6por7eDM0yrrltsholPUQ9/Fmix63NtLqtFgctxiilMWQSsqWxC2Btppdv4KIHh2nTp3Cl19+iT///BPx8fEYMmQIAGDgwIFm7cfOzg6ff/65fnl4U+Tn52PDhg2wtq74ZFhvb2/ExMSgefPmSEhIqPB+LKnSp/a6u7sjI8PwF2pGRgYcHR3L7O9SKpVQKpUWjyVwWBAuHfsGGoVT6YmBEFBoshE4zLJNfFfTrwJoZGI9ehjxmfFGLSL3EhBIL0hHfGY8Orp3rMLIiKgmWbZsGZKSkqBQKODv74/Y2Fj9Mu/m6Nmzp1n1728RmTBhArZs2VJq3aFDh5Y5LdfHxwfvvPOOWceuTJWejHTp0gUHDx40KDty5Ai6dOlS2Yc2olDaQO0nQ9pFAEIYJiT/zMNW+1lBobRwN009GWDCsBdNvcoZGAQUdxXFRW5BQWY2VG7OCJgwFAo74z7M6u56wXWL1iMiul+7du1w8uRJqcMAACxYsAAzZswo9TFHR8cqjqbizE5Gbt26hYsXL+rvJycnIyEhAfXq1UODBg0wZ84cXL16FZ9//jmA4qxt1apVePPNNzFq1Ch8//332LFjBw4cOGC5Z2GGATOfx76le5F+DihSOunLbTQ5cPeTYcDM5y1+zMavBuO3BekosnEus0XGpigLjV8NNn7MAmrTLB5X1YPn3ZtTj4joUebm5ma0imt1ZPaYkV9//RXt2rVDu3btABSv0d+uXTvMmzcPAHDt2jWkpqbq6/v4+ODAgQM4cuQI2rRpgw8++ACfffZZlU/rvdeAmc8jNPJZNO+ohafXTTTvqEVoZP9KSUQAIKBBZ6TU+2cxtftXwvvnfkq9/Qho0Nnix65ts3jau7WHWqWGrIzFkGSQwV3ljvZu7as4MiIiKovZLSMlo23LUtrqqj179sSpU6fMPVSlUiht0Gv001VyLLmVHEGTh+HblevQ7OaLKLpnwTWboiwk1duDPpNHWXwNDM3tO0hObVDuuirJ/8ziqSldNnIrOWYHzEZYTBhkkBkMZC1JUGYFzOJ6I0REj5DaeW0aCQQ3DAYmA0uOL4bPb3VQt8AJWaocJLfOx6wus4sft7C4yC3/ds2U5p91VeIit6D79DEWP75UghsGI6JnBBbHLTYYzKpWqTErYFalnGsiIqo4JiNVKLhhMIK8ghDfs2pWBS3IzLZovepEf665AisR0SOPyUgVk1vJq2xKqcrNGcgysV4NVJXnmoiIKu6RXPSMLCNgwlAoC7OMB82WEALKwpsImPDg5YCJiB5FWp0WJ9JP4OClgziRfqLaLGiYlJQEd3d35OXlmVT/xo0bcHNzw5UrVyo5MmkwGanBFHa28Gnwz8ymMmbx+DRIqzGDV4modnlUr0HVs2dPyGQyo9vdu3f1debMmYPJkyfDwcHBpH26uLhg+PDhCA8PNyuWtLQ0jBo1Cp6enlAoFGjYsCGmTp2Kv//+u8yYbW1t0bRpUyxatKjcCSuWxGSkhus1fzqau5+BUpNtUK7UZKG5+5lKX2dEc/sOji7/DN/NWYajyz/j9XeIyCIe9WtQjR07FteuXTO4lSzhnpqaiv379yM0NNSsfY4cORJbt27FzZs3Tap/6dIldOjQARcuXMCXX36JixcvIjIyEtHR0ejSpYvRfkpiTkpKwpw5czBv3rwyV3C1NCYjtUCv+dMxfE0/tGl2CU3qxqNNs0sYvubZSk9EosOX4/NJB3A6yRcXstrjdJIvPp90oMatbUJEVUvqa1Dl5eVhyJAhqFOnDjw8PLB8+XL07NkT06ZN09dRqVRwd3c3uJXYsWMH2rRpg/r16+vLRo0ahSeeeAKFhYUAAI1Gg3bt2mH48OH6Oi1btoSnpyf27t1rUpyvv/46FAoFvvvuO/To0QMNGjRAnz59EBUVhatXr2Lu3LkG9UtibtiwIUaOHIknnngCR44cqcgpMhuTkVpCYWeL7tPHoPeiGeg+fUyld83UtsXWiKjqmHMNqsoQFhaGY8eOYd++fThy5AhiY2MRH2/6sWJjY9GhQweDso8++gj5+fmYPXs2AGDu3LnIzs7GqlWrDOoFBAQgNjb2gce4efMmDh8+jEmTJhldB87d3R1DhgzB9u3bS+2GEUIgNjYWiYmJUCgUJj+vh8FkhCxOv9gaUPpia4B+sTUiInNJeQ2qvLw8bNq0CcuWLUOvXr3QqlUrbNiwAVqtYSvMmjVrYG9vr7/93//9n/6xlJQUeHp6GtS3t7fHli1bsHr1asybNw8rVqzA5s2bja4v4+npiZSUlAfGeeHCBQgh0KJFi1Ifb9GiBbKysnD9+r/nqCRmpVKJJ598EjqdDlOmTHngsSyBU3vJ4mrrYmtEVDWkvAbVpUuXUFRUhICAAH2Zk5MTmjVrZlBvyJAhBt0gzs7O+v/fvn0btrbGrdNdunTBjBkz8O6772LWrFno3r27UR07OzsUFBSYHO+DBqDe2/JREnNWVhbCw8PRtWtXdO3a1eRjPQwmI2RxtXmxNSKqfCXXoMosyCx13IgMMqhVakmvQeXk5ITGjRuX+piLiwuysowXgdLpdDh27BjkcrnBBWnvdfPmTbi6PjjJaty4MWQyGc6fP4/nnze+7tr58+fh6upqkCTdG/OOHTvQuHFjdO7cGcHBlb9qNbtpyOJMXUStpi62RkSVq+QaVACMLopZ2deg8vX1hY2NDU6cOKEvy8nJwR9//GHyPtq1a4dz584ZlS9duhSJiYn48ccfcejQIWzYsMGoztmzZ/UXqi3PY489hqeffhpr1qzB7du3DR5LT0/H1q1by53NY29vj6lTp2LGjBlVMr2XyQhZHBdbI6LKVnINKjeVm0G5WqVGRM+ISrsGlYODA0aMGIGZM2fihx9+wO+//47Ro0fDysoKsvvHyJUhJCQEx48fNxhncurUKcybNw+fffYZunXrhoiICEydOhWXLl3S1ykoKMDJkyfRu3dvk46zatUqFBYWIiQkBD/99BPS0tJw6NAhPP3002jatCnmzZtX7vbjx4/HH3/8gd27d5t0vIfBZIQsjoutEVFVCG4YjMMvHsb6kPVYErgE60PW49CLhyr9YpgRERHo0qULnn32WQQHB6Nbt25o0aJFqeNAStOnTx9YW1sjKqp4LZQ7d+5g6NChCA0NRf/+/QEA48aNQ1BQEIYNG6ZPWr7++ms0aNAAgYGBJh2nSZMmOHHiBHx9ffHKK6+gYcOG6NOnD5o2bYpjx47B3t6+3O3r1auH4cOH45133oFOpzPpmBUlE1W1vNpDyM3NhZOTE3JycoxGFtOjKzp8OZJTG6BQWVdfpiy8CZ8GaZW+xgkRPdru3LmD5ORk+Pj4mPwl/qjKz89H/fr18cEHH2D06NEmbbN69Wrs27cPhw8fNvk4nTt3xpQpUzB48OCKhorw8HBERETgyJEj6Ny5c4X3c6/yXktTv785gJUqTa/506G5fQdxkVtQkJkNlZszAiYMZYsIEVVrp06dQmJiIgICApCTk4MFCxYAAAYOHGjyPsaPH4/s7Gzk5eWZtCT8jRs38MILL+C1116rcNwAMH/+fHh7e+OXX35BQEAArKwejQ4StowQEVGVq84tI6dOncKYMWOQlJQEhUIBf39/REREoHXr1lUWQ2pqKvz8/Mp8/Ny5c2jQoEGVxMKWEaJyaAqLELslBrnpOXB0d0Lg0J5QKG2kDouIqrl27drh5MmTksbg6emJhISEch+vTpiMUI20b+lepJ8DipROAOrhrzTgz6P74e4HDJhpPOeeiKg6sba2LnMdk+ro0egsIrKgfUv3Iu2iI4oUhk2CRQpHpF10xL6lpl1kioiIqgaTEapRNIVFyDj3zzCoMq6Lk3FOB01hURVHRkREZWEyQjVK7OYfoFE6GyciJWQyaJR1Ebv5hyqNi4iIysZkhGoUzeULFq1HRESVj8kI1SjOjqZ1v5haj4iIKh+TEapR/Ae2gbLw5gOvi+M/sE3VBkZElUJotcj/Xxxy9h9A/v/iIO653kt1MGzYMCxcuNDk+rNnz8bkyZMrMSJpMBmhGkXR5El0UO4svlPGdXE6KHdC0eTJKo6MiCwt97vvcLFXMFJHjMBfM2YgdcQIXOwVjNzvvpMsptatW2PChAmlPrZ582YolUrcuHEDAHD69GkcPHgQU6ZMMXn/M2bMwKZNmwwuoFeey5cvQyaTlbomSUxMDGQyGbKzs40e8/b2xooVK0yO62ExGaGaxUqOtlOGoId2BZSaLIOHlJos9NCuQNspQ4BKuLQ4EVWd3O++w9Wp03A3Pd2g/G5GBq5OnSZZQjJ69Ghs27YNt2/fNnpsw4YNGDBgAFxcXAAAK1euxMsvv/zAC9bdy8XFBSEhIfj4448tFvOjgMkI1Tx+A9AqbCxC/eYhBPPRrWglQjAfoX7z0CpsLOA3QOoIieghCK0WGQsXld4d+09ZxsJFldZlk5eXhyFDhqBOnTrw8PDA8uXL0bNnT0ybNg1Dhw7F7du3sXv3boNtkpOTERMTo7+Qnlarxa5du/RX6QWAxMREqFQqfPHFF/qyHTt2wM7ODufOndOX9e/fH9u2bauU5yYVrsBKNZPfAFg374fGKT8DtzIAezXQsCtbRIhqgIJfTxq1iBgQAnfT01Hw60nU6RRg8eOHhYXh2LFj2LdvH9RqNebNm4f4+Hi0bdsWLi4uGDhwINavX4+hQ4fqt9m4cSMef/xx9O7dGwBw5swZ5OTkoEOHDvo6zZs3x7JlyzBp0iR0794dVlZWmDBhApYsWWJwHZqAgABcuXIFly9fhre3t8WfnxSYjFDNZSUHfAKljoKILOzu9esWrWeOvLw8bNq0CV988QV69eoFoLj75d5rwYwePRp9+vTRXzxOCIFNmzZhxIgR+qvkpqSkQC6Xw83NzWD/kyZNwsGDBzF06FAoFAp07NjRaMBqybFSUlJqTDLCbhoiIqpWrF1dLVrPHJcuXUJRURECAv5tcXFyckKzZs30959++mk8/vjj2LBhAwAgOjoaqampGDlypL7O7du3oVQqIStlgcb169fjzJkziI+Px8aNG43q2NnZAQAKCgos+tykxGSEyMK0Oi1OpJ/AwUsHcSL9BLS66jXVkOhRp+rgD2t393JXWrZ2d4eqg3/VBvYPKysrhIaGYtOmTdDpdNiwYQOCgoLg6+urr+Pi4oKCggJoNBqj7U+fPo38/Hzk5+fj2rVrRo/fvHkTAOD6kMmWo2Px9btycnKMHsvOzoaTk9ND7d8cTEaILCgqJQrP7AjB2pXv4ceV27B25Xt4ZkcIolKipA6NqMaQyeVQvzXnnzulX4NK/dYcyOSWHyPm6+sLGxsbnDhxQl+Wk5ODP/74w6DeyJEjkZaWhj179mDv3r36gasl2rZtCwAGA1OB4kQjNDQUc+fORWhoKIYMGWI0M+fs2bOwsbFBy5YtH+q5NGnSBFZWVjh58qRB+aVLl5CTk4OmTZs+1P7NwTEjRBYSlRKFb1eux8s3J6JIURcA0OA20OG7LHz763pgMhDcMFjiKIlqBsfevYEPVyBj4SKDwazWajXUb80pfrwSODg4YMSIEZg5cybq1asHNzc3hIeHw8rKyqA7xcfHB0899RTGjRsHpVKJF154wWA/rq6uaN++PY4ePapPTABgwoQJ8PLywttvv43CwkK0a9cOM2bMwOrVq/V1YmNjERgYqO+uMUVSUpJRWcuWLTFmzBj83//9H6ytrdG6dWukpaVh1qxZ6Ny5M7p27WrGmXk4TEaILECr0+KHlZvhmzcaRTaGjxXZOMM3bzR+WLkZQe8HQc4ZPUQW4di7Nxx69SqeXXP9OqxdXaHq4F8pLSL3ioiIwIQJE/Dss8/C0dERb775JtLS0mBra2tQb/To0YiOjsakSZOMHgOAMWPG4PPPP8cbb7wBAPj8889x8OBBnDp1CtbW1rC2tsaWLVvQvXt3PPvss+jTpw8AYNu2bXjnnXfMivnVV181KktLS8OHH36IxYsXY9asWUhJSYG7uzuefvpp/Pe//y11PEtlkQlR1rrZj47c3Fw4OTkhJydH38dF9Cg5fvkYfluQjiIb59L7sYWATVEWWs/zQBfvblUeH9Gj5s6dO/rZJqV9UVcn+fn5qF+/Pj744AOj7pjy3L59G82aNcP27dvRpUsXk7b59ttv8X//9384c+YMrK0fjfaE8l5LU7+/OWaEyAIubosq7popZ0BdkaIeLm7j2BGi6u7UqVP48ssv8eeffyI+Ph5DhgwBAAwcONCs/djZ2eHzzz/XLw9vivz8fGzYsOGRSUQspWY9GyKJKG6a1sBoaj0ierQtW7YMSUlJUCgU8Pf3R2xsrH6Zd3P07NnTrPovvfSSwf0JEyZgy5YtpdYdOnQoIiMjzY5JCkxGiCygvnt93Mg1rR4RVW/t2rUzmoEilQULFmDGjBmlPladhjVUqJtm9erV8Pb2hq2tLTp16oS4uLgy6xYVFWHBggVo1KgRbG1t0aZNGxw6dKjCARM9igImDoOiMKv0a2UAgBBQFN5EwMRhVRsYEdVobm5uaNy4cam3+1d3fZSZnYxs374dYWFhCA8PR3x8PNq0aYOQkBBkZmaWWv/tt9/GJ598gpUrV+LcuXOYMGECnn/+eZw6deqhgyd6VCjsbOHbILX4zv0JyT/3fRukQWFXvQfqERFVBrOTkYiICIwdOxYjR46En58fIiMjoVKpsH79+lLrb968GW+99Rb69u0LX19fTJw4EX379sUHH3zw0METPUp6zZ+O5u5noNRkG5QrNVlo7n4GveZPlyYwIqJHnFljRjQaDU6ePIk5c+boy6ysrBAcHIzjx4+Xuk1hYaHRVB87OzscPXq0zOMUFhaisLBQfz8314TOeKJHQK/506G5fQdxkVtQkJkNlZszAiYMZYsIEVE5zEpGbty4Aa1WC7VabVCuVquRmJhY6jYhISGIiIjAk08+iUaNGiE6Ohp79uyBVlv29ToWLVqE+fPnmxMa0SNDYWeL7tPHSB0GEVG1UenrjHz44Ydo0qQJmjdvDoVCgTfeeAMjR47UX0a5NHPmzEFOTo7+lpaWVtlhEhERkUTMSkZcXFwgl8uRkZFhUJ6RkQF3d/dSt3F1dcVXX32F/Px8pKSkIDExEfb29gZXL7yfUqmEo6OjwY2IiKimSEpKgru7O/Ly8kyqf+PGDbi5ueHKlSuVHJk0zEpGShZ3iY6O1pfpdDpER0c/cClbW1tb1K9fH3fv3sXu3bvNXqmOiIjofjqdwNWkLPxxIh1Xk7Kg00m/sGDPnj0hk8mMbnfv3tXXmTNnDiZPngwHBweT9uni4oLhw4cjPDzc5DhCQ0Px3HPPlRnjtGnTjMo3btwIZ2dnk49hKWYvehYWFoYRI0agQ4cOCAgIwIoVK5Cfn4+RI0cCAIYPH4769etj0aJFAID//e9/uHr1Ktq2bYurV6/inXfegU6nw5tvvmnZZ0JERLXKn6cyEbv9AvKz/53wUMdZicBBTdConbRrbIwdOxYLFiwwKCtZwj01NRX79+/HypUrzdrnyJEj4e/vj6VLl6JevXoWi/VRYPaYkUGDBmHZsmWYN28e2rZti4SEBBw6dEg/qDU1NRXXrl3T179z5w7efvtt+Pn54fnnn0f9+vVx9OhRSTIvIiKqGf48lYlDn5w1SEQAID+7EIc+OYs/T5W+9pUl5OXlYciQIahTpw48PDywfPlyo5YGlUoFd3d3g1uJHTt2oE2bNqhf/98VmUeNGoUnnnhCP5NUo9GgXbt2GD58uL5Oy5Yt4enpib1791bac5NKhZaDf+ONN/SXPL5fTEyMwf0ePXrg3LlzFTkMEZmBU4qpttDpBGK3Xyi3ztEdF+DTxhVWVmVcvPIhhIWF4dixY9i3bx/UajXmzZuH+Ph4tG3b1qTtY2Nj0aFDB4Oyjz76CG3atMHs2bOxfPlyzJ07F9nZ2Vi1apVBvYCAAMTGxpp1deDqgNemIaoBosOXIzm1AQqV/wwMzwISJx2AT4NULrZGNc61C9lGLSL3u5VViGsXslG/WV2LHjsvLw+bNm3CF198gV69egEANmzYAE9PT4N6a9aswWeffaa/P378eP1inykpKUbJiL29PbZs2YIePXrAwcEBK1aswA8//GA0gcPT07NGrmDOZISomosOX47E9CcAhWF5ocIZienOQPhyJiRUo+Tnlp+ImFvPHJcuXUJRURECAgL0ZU5OTmjWrJlBvSFDhmDu3Ln6+/cOTbh9+7bRYqAA0KVLF8yYMQPvvvsuZs2ahe7duxvVsbOzQ0FBgQWeyaOFyQhRNaa5fQfJqQ2KExHZfc3RMhkgBJJTvaC5fYddNlRj1HFUWrReZXByckLjxo1LfczFxQVZWVlG5TqdDseOHYNcLsfFixdL3fbmzZtwdXV96PgcHR2Rk5NjVJ6dnQ0nJ6eH3r+5Kn3RMyKqPHGRW1CorGuciJSQyVCorIe4yC1VGxhRJfJo4ow6zuUnGvZ1lfBo4mzxY/v6+sLGxgYnTpzQl+Xk5OCPP/4weR/t2rUrdSzl0qVLkZiYiB9//BGHDh3Chg0bjOqcPXsW7dq1q1jw92jWrBni4+ONyuPj49G0adOH3r+5mIwQVWMFmdkWrUdUHVhZyRA4qEm5dbq/0qRSBq86ODhgxIgRmDlzJn744Qf8/vvvGD16NKysrCAr60fBfUJCQnD8+HGDy6KcOnUK8+bNw2effYZu3bohIiICU6dOxaVLl/R1CgoKcPLkSfTu3dvkeHNycpCQkGBwS0tLw8SJE/HHH39gypQpOHPmDJKSkhAREYEvv/wS//d//2f6CbEQJiNE1ZjKzdmi9Yiqi0bt3PDM+FZGLST2dZV4ZnyrSl1nJCIiAl26dMGzzz6L4OBgdOvWDS1atCh1HEhp+vTpA2tra0RFRQEoXgJj6NChCA0NRf/+/QEA48aNQ1BQEIYNG6ZPWr7++ms0aNAAgYGBJscaExODdu3aGdzmz58PX19f/PTTT0hMTERwcDA6deqEHTt2YOfOnXjmmWfMPCMPTyaEkH65ugfIzc2Fk5MTcnJyuDQ80T00t+/g80kHUKhwLr2rRggoNVkYvuZZjhmhR8qdO3eQnJwMHx8fk7/ES6PTieLZNbmFqONY3DVTGS0i5cnPz0f9+vXxwQcfmDzldvXq1di3bx8OHz5s8nE6d+6MKVOmYPDgwRUNtVKU91qa+v3NAaxE1ZjCzhY+DVKLZ80IYZiQ/PM7w6dBGhMRqrGsrGQWn777IKdOnUJiYiICAgKQk5OjX2nVnMucjB8/HtnZ2cjLyzNpSfgbN27ghRdewGuvvVbhuB9l7KYhquZ6zZ+O5u5noNRkG5QrNVlo7n6G03qJKsGyZcvQpk0bBAcHIz8/H7GxsXBxcTF5e2tra8ydO9esa9O8+eab+nEpqampsLe3L/OWmppaoeclFbaMENUAveZP5wqsRFWkXbt2OHnypKQxeHp6IiEhodzHqxMmI0Q1hMLOFt2nj5E6DCKqAtbW1mWuY1IdsZuGiIgkUw3mUNADWOI1ZDJCRERVzsbGBgBq5NLmtU3Ja1jymlYEu2mIiKjKyeVyODs7IzMzEwCgUqlMXjSMHg1CCBQUFCAzMxPOzs6Qy+UV3heTESIikoS7uzsA6BMSqp6cnZ31r2VF1dpkRKvTIj4zHtcLrsNV5Yr2bu0ht6p4VkdEROaRyWTw8PCAm5sbioqKpA6HKsDGxuahWkRK1MpkJColCovjFiOjIENfplapMTtgNoIbBksYGRFR7SOXyy3yhUbVV60bwBqVEoWwmDCDRAQAMgsyERYThqiUKIkiI6qmdFogORb4bVfxvzrtg7chIrpHrWoZ0eq0WBy3GALG05AEBGSQYUncEgR5BbHLhsgU5/bh7oHZuPynK24V1YO9zU14N7oO636LAb8BUkdHRNVErUpG4jPjjVpE7iUgkF6QjvjMeHR071iFkRFVQ+f24WzEWvyiXYBCZT3gn1l9ynM30TlpLVqFgQkJEZmkVnXTXC+4btF6RLWWTouEj7biR/k0FCoML1JWqKiLH+XTkPDRVnbZEJFJalUy4qpytWg9otpKc+En/Fr4cvGd+9eG+Of+r4UvQ3PhpyqOjIiqo1qVjLR3aw+1Sg0ZSl9YRwYZ3FXuaO/WvoojI3p4Wp0WJ9JP4OClgziRfgLaSmyVOPn16eKumbIWqZLJUKish5Nfn660GIio5qhVY0bkVnLMDpiNsJgwyCAzGMhakqDMCpjFwatU7VT1dPXsXNOWfTa1HhHVbrWqZQQAghsGI6JnBNxUbgblapUaET0juM4IVTtSTFdXeDexaD0iqt1kohpcMjE3NxdOTk7IycmBo6OjRfbJFVipJtDqtAjZHVLmLDEZZFCr1Dj04iGLvr81hUXYNOEbaBROpXfVCAGFJhsjIgdAoWTrCFFtZer3d63qprmX3ErO6btU7Uk1XV2htIHaT4a0iwCEMExI/vl9o/azYiJCRCapdd00RDWJlNPVB8x8Hl6Nc2GjyTUot9HkwKtxLgbMfN7ixySimqnWtowQ1QRST1cfMPN5aAqLELslBrnpOXB0d0Lg0P5sESEiszAZIarGSqarZxZklnqZg5IxI5U5XV2htEGv0U9X2v6JqOZjNw1RNVYyXR2A0fo5nK5ORNUFkxGiao7T1YmoumM3DVENENwwGEFeQZyuTkTVEpMRohqC09WJqLpiNw0RERFJiskIERERSYrJCBEREUmKY0aIqFrS3L6DuMgtKMjMhsrNGQEThkJhZyt1WERUARVqGVm9ejW8vb1ha2uLTp06IS4urtz6K1asQLNmzWBnZwcvLy9Mnz4dd+7cqVDARETR4cvx+aQDOJ3kiwtZ7XE6yRefTzqA6PDlUodGRBVgdjKyfft2hIWFITw8HPHx8WjTpg1CQkKQmZlZav0vvvgCs2fPRnh4OM6fP49169Zh+/bteOuttx46eCKqfaLDlyMx/QkUKpwNygsVzkhMf4IJCVE1ZHYyEhERgbFjx2LkyJHw8/NDZGQkVCoV1q9fX2r9n3/+Gd26dcPgwYPh7e2N3r1747XXXntgawoR0f00t+8gObVB8R2Z4YqzJfeTU72guc2WV6LqxKxkRKPR4OTJkwgO/ndFRysrKwQHB+P48eOlbtO1a1ecPHlSn3xcunQJBw8eRN++fcs8TmFhIXJzcw1uRERxkVtQqKxrnIiUkMlQqKyHuMgtVRsYET0Uswaw3rhxA1qtFmq12qBcrVYjMTGx1G0GDx6MGzduoHv37hBC4O7du5gwYUK53TSLFi3C/PnzzQmNiGqBgsxsi9YjokdDpU/tjYmJwcKFC7FmzRrEx8djz549OHDgAN59990yt5kzZw5ycnL0t7S0tMoOk4iqAZWbs0XrEdGjwayWERcXF8jlcmRkZBiUZ2RkwN3dvdRt/vOf/2DYsGEYM2YMAKB169bIz8/HuHHjMHfuXFhZGedDSqUSSqXSnNCIqBYImDAUiZMOFA9eLa2rRggoNVkImDC0ymMjooozq2VEoVDA398f0dHR+jKdTofo6Gh06dKl1G0KCgqMEg65vPjiXUIIc+MlolpMYWcLnwapxXfu//vxz32fBmlcb4SomjG7myYsLAxr167Fpk2bcP78eUycOBH5+fkYOXIkAGD48OGYM2eOvn7//v3x8ccfY9u2bUhOTsaRI0fwn//8B/3799cnJUREpuo1fzqau5+BUpNtUK7UZKG5+xn0mj9dmsCIqMLMXoF10KBBuH79OubNm4f09HS0bdsWhw4d0g9qTU1NNWgJefvttyGTyfD222/j6tWrcHV1Rf/+/fHf//7Xcs+CiGqVXvOncwVWohpEJqpBX0lubi6cnJyQk5MDR0dHqcMhIiIiE5j6/c0L5REREZGkmIwQERGRpJiMEBERkaSYjBAREZGkmIwQERGRpMye2ktEdC+tTov4zHhcL7gOV5Ur2ru1h9yKawgRkemYjBBRhUWlRGFx3GJkFPx7iQi1So3ZAbMR3DC4nC2JiP7FbhoiqpColCiExYQZJCIAkFmQibCYMESlREkUGRFVN2wZISKzaXVaLI5bDAHjNRMFBGSQYUncEgR5BdW4Lhuu/EpkeUxGiMhs8ZnxRi0i9xIQSC9IR3xmPDq6d6zCyCpXdPhyJKc2QKHSt7ggC0icdAA+DVJ5TRyih8BuGiIy2/WC6xatVx1Ehy9HYvoTKFQ4G5QXKpyRmP4EosOXSxMYUQ3AZISIzOaqcrVovUed5vYdJKc2KL4jkxk++M/95FQvaG7fqeLIiGoGJiNEZLb2bu2hVqkhg6zUx2WQwV3ljvZu7as4ssoRF7kFhcq6xolICZkMhcp6iIvcUrWBEdUQTEaIyGxyKzlmB8wGAKOEpOT+rIBZNWbwakFmtkXrEZEhJiNEVCHBDYMR0TMCbio3g3K1So2InhE1ap0RlZuzResRkSGZEMJ4bt4jJjc3F05OTsjJyYGjo6PU4RDRPWrDCqya23fw+aQDxYNXS+uqEQJKTRaGr3mW03yJ7mHq9zen9hLRQ5FbyWvU9N3SKOxs4dMgFYnpzoAQhgnJP7/nfBqkMREhqiB20xARmaDX/Olo7n4GSk22QblSk4Xm7me4zgjRQ2DLCBGRiXrNn84VWIkqAZMRIiIzKOxs0X36GKnDIKpR2E1DREREkmIyQkRERJJiNw0REZWpNkzdJukxGSEiolJFpURhcdxigys0q1VqzA6YXaMWtSPpsZuGiIiMRKVEISwmzCARAYDMgkyExYQhKiVKosioJmIyQkTVklanxYn0Ezh46SBOpJ+AVqeVOqQaQ6vTYnHcYggYL9BdUrYkbgnPOVkMu2mIqNqpjd0HVTl2Iz4z3qhF5F4CAukF6YjPjK+01Xc5VqV2YTJCRNVKSffB/b/aS7oPatpF+oCqT76uF1y3aD1z1cZks7ZjNw0RVRu1sftAirEbripXi9YzB8eq1E5MRoio2jCn+6AmuDf5stICXf7wRd+Edujyhy9kWgFAVEry1d6tPdQqNWQo5QrFAGSQwV3ljvZu7S163NqYbFIxdtMQUbUhdfcBIM3Yjd6/tUazmy+iSFEXANDgNtDhWhaS6u3Gd61/s/jYDbmVHLMDZiMsJgwyyAySg5IEZVbALIs/70dhrApJg8kIEVUbUnYfABKM3cgvTkR880ajyMbwsSIbZ/jmjUbv39bheveyv8ArKrhhMCJ6RpT6fGcFzKqRY1VIOkxGiKjaKOk+yCzILLUpXwYZ1Cq1xbsPAGkGztb9+6/iFhEbALL7ukxkMkAINLv5Aur+/RfQyKKHBlCckAR5BVVZS5DUySZJh2NGiKjaKOk+AGA0nqEyuw+kGssgYm4Vd83cn4iUkMlQpKgHEXPLose9l9xKjo7uHdHXty86unes1Om1Uo1VIekxGSGiaqWk+8BN5WZQrlapK21ar1QDZ3NylRat96iTKtkk6bGbhoiqnaruPpBqLIPCuwlww8R6NYQUY1VIekxGiKhaKuk+qApSjWUIHBaES8e+gUbhVHpXjRBQaLIROGyARY8rtapONkl6TEaIiB5AqoGzCqUN1H4ypF0EIIRhQiKK41D7WUGhtCl9B9VYVSabJD2OGSEiegApxzIMmPk8vBrnwkaTa1Buo8mBV+NcDJj5vMWPSVTVKpSMrF69Gt7e3rC1tUWnTp0QFxdXZt2ePXtCJpMZ3fr161fhoImIqpoUA2dLDJj5PEIjn0Xzjlp4et1E845ahEb2ZyJCNYZMCGHc5liO7du3Y/jw4YiMjESnTp2wYsUK7Ny5E0lJSXBzczOqf/PmTWg0Gv39v//+G23atMFnn32G0NBQk46Zm5sLJycn5OTkwNHR0ZxwiYgsileTJTKdqd/fZicjnTp1QseOHbFq1SoAgE6ng5eXFyZPnozZs2c/cPsVK1Zg3rx5uHbtGurUqWPSMZmMEBERVT+mfn+b1U2j0Whw8uRJBAf/2xxpZWWF4OBgHD9+3KR9rFu3Dq+++mq5iUhhYSFyc3MNbkRERFQzmZWM3LhxA1qtFmq12qBcrVYjPT39gdvHxcXh7NmzGDNmTLn1Fi1aBCcnJ/3Ny8vLnDCJiIioGqnS2TTr1q1D69atERAQUG69OXPmICcnR39LS0urogiJiIioqpm1zoiLiwvkcjkyMgyXRc7IyIC7u3u52+bn52Pbtm1YsGDBA4+jVCqhVNaM5Y2JiIiofGa1jCgUCvj7+yM6OlpfptPpEB0djS5dupS77c6dO1FYWIihQ4dWLFIiIiKqkcxegTUsLAwjRoxAhw4dEBAQgBUrViA/Px8jR44EAAwfPhz169fHokWLDLZbt24dnnvuOTz22GOWiZyIiIhqBLOTkUGDBuH69euYN28e0tPT0bZtWxw6dEg/qDU1NRVWVoYNLklJSTh69Ci+++47y0RNRERVQlNYhNgtMchNz4GjuxMCh/askcvPk7TMXmdEClxnhIio6u1buhfp54AipZO+zKYwB+5+4OqvZJJKWWeEiIhqh31L9yLtoiOKFIZfIEUKR6RddMS+pXslioxqIiYjRERkQFNYhIxz/zSaywwvDFhyP+OcDprCoiqOjGoqJiNERGQgdvMP0CidjROREjIZNMq6iN38Q5XGVRW0Oi1OpJ/AwUsHcSL9BLQ6rdQh1QpmD2AlIqKaTXP5AoAWJtbrXenxVJWolCgsjluMjIJ/19JSq9SYHTC7Uq/KTGwZISKi+zg7mtb9Ymq96iAqJQphMWEGiQgAZBZkIiwmDFEpURJFVjswGSEiIgP+A9tAWXgTKGuypRBQFt6E/8A2VRtYJdHqtFgctxgCxs+3pGxJ3BJ22VQiJiNERGRA0eRJdFDuLL5zf0Lyz/0Oyp1QNHmyiiOrHPGZ8UYtIvcSEEgvSEd8ZnwVRlW7MBkhIiJDVnK0nTIEPbQroNRkGTyk1GShh3YF2k4ZAljJJQrQsq4XXLdoPTIfB7ASEZExvwFoFQY0PzAbl/90xa2ierC3uQlvv+uw7rcY8BsgdYQW46pytWg9Mh+TESIiKp3fAFg374fGKT8DtzIAezXQsGuNaREp0d6tPdQqNTILMksdNyKDDGqVGu3d2ksQXe3AbhoiIiqblRzwCQRav1T8bw1LRABAbiXH7IDZAIoTj3uV3J8VMAvyGvjcHxVMRoiIqNYLbhiMiJ4RcFO5GZSrVWpE9IzgOiOVjN00REREKE5IgryCEJ8Zj+sF1+GqckV7t/ZsEakCTEaIiIj+IbeSo6N7R6nDqHXYTUNERESSYjJCREREkmIyQkRERJJiMkJERESS4gBWIiJ65Ghu30Fc5BYUZGZD5eaMgAlDobCzlTosqiRMRoiI6JESHb4cyakNUKj0LS7IAhInHYBPg1T0mj9d2uCoUrCbhoiIHhnR4cuRmP4EChXOBuWFCmckpj+B6PDl0gRGlYrJCBERPRI0t+8gObVB8R2Z4bLsJfeTU72guX2niiOjysZkhIiIHglxkVtQqKxrnIiUkMlQqKyHuMgtVRsYVTomI0RE9EgoyMy2aD2qPpiMEBHRI0Hl5mzRelR9MBkhIqJHQsCEoVAWZgFClF5BCCgLbyJgwtCqDYwqHZMRIiJ6JCjsbOHTILX4zv0JyT/3fRqkcb2RGojJCBERPTJ6zZ+O5u5noNRkG5QrNVlo7n6G64zUUFz0jIiIHim95k/nCqy1DJMRIiJ65CjsbNF9+hipw6Aqwm4aIiIikhSTESIiIpIUkxEiIiKSFJMRIiIikhSTESIiIpIUkxEiIiKSFJMRIiIikhTXGSEiIvoHF1uTBpMRIiIiANHhy5Gc2gCFSt/igiwgcdIB+DRI5TL0laxC3TSrV6+Gt7c3bG1t0alTJ8TFxZVbPzs7G6+//jo8PDygVCrRtGlTHDx4sEIBExERWVp0+HIkpj+BQoWzQXmhwhmJ6U8gOny5NIHVEma3jGzfvh1hYWGIjIxEp06dsGLFCoSEhCApKQlubm5G9TUaDZ5++mm4ublh165dqF+/PlJSUuDs7GyJ+ImIiB6K5vYdJKc2ABQAZDLDB2UyQAgkp3pBc/sOu2wqiUyI+6/TXL5OnTqhY8eOWLVqFQBAp9PBy8sLkydPxuzZs43qR0ZGYunSpUhMTISNjU2FgszNzYWTkxNycnLg6OhYoX0QERGV5ujyz3A6yfeB9do0u8Tr5ZjJ1O9vs7ppNBoNTp48ieDg4H93YGWF4OBgHD9+vNRt9u3bhy5duuD111+HWq1Gq1atsHDhQmi12jKPU1hYiNzcXIMbERFRZSjIzLZoPTKfWcnIjRs3oNVqoVarDcrVajXS09NL3ebSpUvYtWsXtFotDh48iP/85z/44IMP8N5775V5nEWLFsHJyUl/8/LyMidMIiIik6ncnC1aj8xX6euM6HQ6uLm54dNPP4W/vz8GDRqEuXPnIjIyssxt5syZg5ycHP0tLS2tssMkIqJaKmDCUCgLs4CyRi0IAWXhTQRMGFq1gdUiZg1gdXFxgVwuR0ZGhkF5RkYG3N3dS93Gw8MDNjY2kMvl+rIWLVogPT0dGo0GCoXCaBulUgmlUmlOaERERBWisLOFT4NUJKY7Fyck9w5i/SdB8WmQxsGrlcislhGFQgF/f39ER0fry3Q6HaKjo9GlS5dSt+nWrRsuXrwInU6nL/vjjz/g4eFRaiJCRERU1XrNn47m7meg1GQblCs1WWjufobrjFQys6f2hoWFYcSIEejQoQMCAgKwYsUK5OfnY+TIkQCA4cOHo379+li0aBEAYOLEiVi1ahWmTp2KyZMn48KFC1i4cCGmTJli2WdCRET0EHrNn84VWCVidjIyaNAgXL9+HfPmzUN6ejratm2LQ4cO6Qe1pqamwsrq3wYXLy8vHD58GNOnT8cTTzyB+vXrY+rUqZg1a5blngUREZEFKOxsOX1XAmavMyIFrjNCRERU/VTKOiNERERElsZkhIiIiCTFZISIiIgkxWSEiIiIJGX2bBoiIiKyLE1hEWK3xCA3PQeO7k4IHNoTCmXFLi5bHTEZISIiktC+pXuRfg4oUjoBqIe/0oA/j+6Hux8wYObzUodXJdhNQ0REJJF9S/ci7aIjihSG016LFI5Iu+iIfUv3ShRZ1WIyQkREJAFNYREyzv2z1Ne918O5537GOR00hUVVHFnVYzJCREQkgdjNP0CjdDZORErIZNAo6yJ28w9VGpcUmIwQERFJQHP5gkXrVWdMRoiIiCTg5Fho0XrVGZMRIiIiCch62sNGkwWUdYk4IWCjuQlZT/uqDEsSTEaIiIgkkPWYJ5Lq7S6+c39C8s/9pHp7kPWYZxVHVvWYjBAREUnAtY4a37X+DZcc1sGmKNvgMZuiLFxyWIfvWv8G1zpqaQKsQlz0jIiISALt3dpDrVLjSOvfEKX9DZ3+9EXdAidkqXLwv0aXIOSAu8od7d3aSx1qpWMyQkREJAG5lRyzA2YjLCYMQg4cb3pJ/5gMxdN9ZwXMgtxKLlWIVYbdNERERBIJbhiMiJ4RcFO5GZSrVWpE9IxAcMNgiSKrWmwZISIiklBww2AEeQUhPjMe1wuuw1XlivZu7WtFi0gJJiNEREQSk1vJ0dG9o9RhSIbdNERERCQptowQERHVVjotkPIzcCsDsFcDDbsCEnQPMRkhIiKqjc7tw90Ds3H5T1fcKqoHe5ub8G50Hdb9FgN+A6o0FCYjREREtc25fTgbsRa/aBegUFkPsCkuVp67ic5Ja9EqDFWakHDMCBERUW2i0yLho634UT4NhYq6Bg8VKuriR/k0JHy0tbgLp4owGSEiIqpFNBd+wq+FLxffkckMH/zn/q+FL0Nz4acqi4nJCBERUS1y8uvTxV0z9yciJWQyFCrr4eTXp6ssJiYjREREtUh2ro1F61kCkxEiIqJaROHdxKL1LIHJCBERUS0SOCwIisJsQIjSKwgBRWEWAocFVVlMTEaIiIhqEYXSBmq/f8aL3J+Q/HNf7WcFhZLdNERERFRJBsx8Hl6Nc2GjyTUot9HkwKtxLgbMfL5K4+GiZ0RERLXQgJnPQ1NYhNgtMchNz4GjuxMCh/av0haREkxGiIiIaimF0ga9Rj8tdRjspiEiIiJpMRkhIiIiSTEZISIiIkkxGSEiIiJJMRkhIiIiSTEZISIiIklVKBlZvXo1vL29YWtri06dOiEuLq7Muhs3boRMJjO42draVjhgIiIiqlnMTka2b9+OsLAwhIeHIz4+Hm3atEFISAgyMzPL3MbR0RHXrl3T31JSUh4qaCIiIqo5zE5GIiIiMHbsWIwcORJ+fn6IjIyESqXC+vXry9xGJpPB3d1df1Or1eUeo7CwELm5uQY3IiIiqpnMSkY0Gg1OnjyJ4ODgf3dgZYXg4GAcP368zO1u3bqFhg0bwsvLCwMHDsTvv/9e7nEWLVoEJycn/c3Ly8ucMImIiKgaMSsZuXHjBrRarVHLhlqtRnp6eqnbNGvWDOvXr8fXX3+NLVu2QKfToWvXrrhy5UqZx5kzZw5ycnL0t7S0NHPCJCIiomqk0q9N06VLF3Tp0kV/v2vXrmjRogU++eQTvPvuu6Vuo1QqoVQqKzs0IiIiegSY1TLi4uICuVyOjIwMg/KMjAy4u7ubtA8bGxu0a9cOFy9eNOfQREREVEOZlYwoFAr4+/sjOjpaX6bT6RAdHW3Q+lEerVaL3377DR4eHuZFSkRERDWS2d00YWFhGDFiBDp06ICAgACsWLEC+fn5GDlyJABg+PDhqF+/PhYtWgQAWLBgATp37ozGjRsjOzsbS5cuRUpKCsaMGWPZZ0JERETVktnJyKBBg3D9+nXMmzcP6enpaNu2LQ4dOqQf1Jqamgorq38bXLKysjB27Fikp6ejbt268Pf3x88//ww/Pz/LPQsiIiKqtmRCCCF1EA+Sm5sLJycn5OTkwNHRUepwiIiIyASmfn/z2jREREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJKkKJSOrV6+Gt7c3bG1t0alTJ8TFxZm03bZt2yCTyfDcc89V5LBERERUA5mdjGzfvh1hYWEIDw9HfHw82rRpg5CQEGRmZpa73eXLlzFjxgwEBgZWOFgiIiKqecxORiIiIjB27FiMHDkSfn5+iIyMhEqlwvr168vcRqvVYsiQIZg/fz58fX0fKmAiIiKqWcxKRjQaDU6ePIng4OB/d2BlheDgYBw/frzM7RYsWAA3NzeMHj3apOMUFhYiNzfX4EZEREQ1k1nJyI0bN6DVaqFWqw3K1Wo10tPTS93m6NGjWLduHdauXWvycRYtWgQnJyf9zcvLy5wwiYiIqBqp1Nk0eXl5GDZsGNauXQsXFxeTt5szZw5ycnL0t7S0tEqMkoiIiKRkbU5lFxcXyOVyZGRkGJRnZGTA3d3dqP6ff/6Jy5cvo3///voynU5XfGBrayQlJaFRo0ZG2ymVSiiVSnNCIyIiomrKrJYRhUIBf39/REdH68t0Oh2io6PRpUsXo/rNmzfHb7/9hoSEBP1twIABCAoKQkJCArtfiIiIyLyWEQAICwvDiBEj0KFDBwQEBGDFihXIz8/HyJEjAQDDhw9H/fr1sWjRItja2qJVq1YG2zs7OwOAUTkRERFVLa1Oi/jMeFwvuA5XlSvau7WH3Epe5XGYnYwMGjQI169fx7x585Ceno62bdvi0KFD+kGtqampsLLiwq5ERESPsqiUKCyOW4yMgn+HXqhVaswOmI3ghsHlbGl5MiGEqNIjVkBubi6cnJyQk5MDR0dHqcMhIiKq1qJSohAWEwYBwxRABhkAIKJnhEUSElO/v9mEQUREVItodVosjltslIgA0JctiVsCrU5bZTExGSEiIqpF4jPjDbpm7icgkF6QjvjM+CqLickIERFRLXK94LpF61kCkxEiIqJaxFXlatF6lsBkhIiIqBZp79YeapVaP1j1fjLI4K5yR3u39lUWE5MRIiKiWkRuJcfsgNkAYJSQlNyfFTCrStcbYTJCRERUywQ3DEZEzwi4qdwMytUqtcWm9ZrD7EXPiIiIqPoLbhiMIK+g6rkCKxEREdUMcis5Orp3lDoMdtMQERGRtJiMEBERkaSYjBAREZGkmIwQERGRpJiMEBERkaSYjBAREZGkmIwQERGRpJiMEBERkaSYjBAREZGkqsUKrEIIAEBubq7EkRAREZGpSr63S77Hy1ItkpG8vDwAgJeXl8SREBERkbny8vLg5ORU5uMy8aB05RGg0+nw119/wcHBATKZ7MEbmCg3NxdeXl5IS0uDo6OjxfZLpuH5lx5fA2nx/EuL57/yCSGQl5cHT09PWFmVPTKkWrSMWFlZ4fHHH6+0/Ts6OvKNKCGef+nxNZAWz7+0eP4rV3ktIiU4gJWIiIgkxWSEiIiIJFWrkxGlUonw8HAolUqpQ6mVeP6lx9dAWjz/0uL5f3RUiwGsREREVHPV6pYRIiIikh6TESIiIpIUkxEiIiKSFJMRIiIikhSTESIiIpJUjU9GVq9eDW9vb9ja2qJTp06Ii4srt/7OnTvRvHlz2NraonXr1jh48GAVRVozmXP+165di8DAQNStWxd169ZFcHDwA18vejBzPwMltm3bBplMhueee65yA6zhzD3/2dnZeP311+Hh4QGlUommTZvy79BDMPf8r1ixAs2aNYOdnR28vLwwffp03Llzp4qircVEDbZt2zahUCjE+vXrxe+//y7Gjh0rnJ2dRUZGRqn1jx07JuRyuXj//ffFuXPnxNtvvy1sbGzEb7/9VsWR1wzmnv/BgweL1atXi1OnTonz58+L0NBQ4eTkJK5cuVLFkdcc5r4GJZKTk0X9+vVFYGCgGDhwYNUEWwOZe/4LCwtFhw4dRN++fcXRo0dFcnKyiImJEQkJCVUcec1g7vnfunWrUCqVYuvWrSI5OVkcPnxYeHh4iOnTp1dx5LVPjU5GAgICxOuvv66/r9Vqhaenp1i0aFGp9V955RXRr18/g7JOnTqJ8ePHV2qcNZW55/9+d+/eFQ4ODmLTpk2VFWKNV5HX4O7du6Jr167is88+EyNGjGAy8hDMPf8ff/yx8PX1FRqNpqpCrNHMPf+vv/66eOqppwzKwsLCRLdu3So1ThKixnbTaDQanDx5EsHBwfoyKysrBAcH4/jx46Vuc/z4cYP6ABASElJmfSpbRc7//QoKClBUVIR69epVVpg1WkVfgwULFsDNzQ2jR4+uijBrrIqc/3379qFLly54/fXXoVar0apVKyxcuBBarbaqwq4xKnL+u3btipMnT+q7ci5duoSDBw+ib9++VRJzbVYtrtpbETdu3IBWq4VarTYoV6vVSExMLHWb9PT0Uuunp6dXWpw1VUXO//1mzZoFT09PowSRTFOR1+Do0aNYt24dEhISqiDCmq0i5//SpUv4/vvvMWTIEBw8eBAXL17EpEmTUFRUhPDw8KoIu8aoyPkfPHgwbty4ge7du0MIgbt372LChAl46623qiLkWq3GtoxQ9bZ48WJs27YNe/fuha2trdTh1Ap5eXkYNmwY1q5dCxcXF6nDqZV0Oh3c3Nzw6aefwt/fH4MGDcLcuXMRGRkpdWi1QkxMDBYuXIg1a9YgPj4ee/bswYEDB/Duu+9KHVqNV2NbRlxcXCCXy5GRkWFQnpGRAXd391K3cXd3N6s+la0i57/EsmXLsHjxYkRFReGJJ56ozDBrNHNfgz///BOXL19G//799WU6nQ4AYG1tjaSkJDRq1Khyg65BKvIZ8PDwgI2NDeRyub6sRYsWSE9Ph0ajgUKhqNSYa5KKnP///Oc/GDZsGMaMGQMAaN26NfLz8zFu3DjMnTsXVlb8/V5ZauyZVSgU8Pf3R3R0tL5Mp9MhOjoaXbp0KXWbLl26GNQHgCNHjpRZn8pWkfMPAO+//z7effddHDp0CB06dKiKUGssc1+D5s2b47fffkNCQoL+NmDAAAQFBSEhIQFeXl5VGX61V5HPQLdu3XDx4kV9EggAf/zxBzw8PJiImKki57+goMAo4ShJDAWvKVu5pB5BW5m2bdsmlEql2Lhxozh37pwYN26ccHZ2Funp6UIIIYYNGyZmz56tr3/s2DFhbW0tli1bJs6fPy/Cw8M5tfchmHv+Fy9eLBQKhdi1a5e4du2a/paXlyfVU6j2zH0N7sfZNA/H3POfmpoqHBwcxBtvvCGSkpLE/v37hZubm3jvvfekegrVmrnnPzw8XDg4OIgvv/xSXLp0SXz33XeiUaNG4pVXXpHqKdQaNToZEUKIlStXigYNGgiFQiECAgLEL7/8on+sR48eYsSIEQb1d+zYIZo2bSoUCoVo2bKlOHDgQBVHXLOYc/4bNmwoABjdwsPDqz7wGsTcz8C9mIw8PHPP/88//yw6deoklEql8PX1Ff/973/F3bt3qzjqmsOc819UVCTeeecd0ahRI2Frayu8vLzEpEmTRFZWVtUHXsvIhGDbExEREUmnxo4ZISIiouqByQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJ6v8B4qMWGZVUQAEAAAAASUVORK5CYII="
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGzCAYAAAD9pBdvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABnv0lEQVR4nO3deVxUVf8H8M8wMMMiq+wqAu77PkpqapKUpvbYk5ai4L5rkqWmSdqTSylaLlnmlvtambtQFi5FiWimoAWCG4uyBsiwnN8fyPwa2WZw4AJ+3q/XvPSee+6937mz8J1zzj1XJoQQICIiIpKIkdQBEBER0bONyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCNUaH3zwAWQymV51Hzx4UMlR1UxnzpyBTCbDmTNnpA6lTP7+/nB3d5c6jBrL3d0d/v7+OtXt3bs3evfuXanx0LOrRiUjW7duhUwmw++//y51KFquXbuGDz74ALdu3dJ5m7Nnz+Lll19GvXr1YGpqCjc3NwwcOBC7du2qvECfQUuWLMG3335bafs/cuQIXnrpJdStWxempqZo2rQp3nnnHSQnJxer6+/vD5lMpnkolUo0bdoUCxcuxKNHjyotxspU3mfy1q1bkMlkWLFiRYnrV6xYAZlMptdnhypPRb7LKltRYiyTybBjx44S63Tv3h0ymQytW7eu0tju37+PuXPnok+fPrC0tCw3gVer1ViyZAmaN28OU1NTODk5YcCAAbhz506Zxyn6HMlkMhw8eLDYeql+XN2/fx8TJkyAh4cHzMzM0KhRIwQEBODhw4d676tGJSPV1bVr17Bo0SKdP8D79+/H888/j4SEBMycORNr1qyBr68vUlJSsHHjxsoNthZbsGABsrOztcoqMxmZPXs2Bg4ciPj4eMyZMwdr166Ft7c31qxZg/bt2+PmzZvFtlEqldi+fTu2b9+OoKAguLu748MPP8TYsWMrJUaiskRFRWl955T1XXbq1CmcOnWqCqPTZmpqWuKPtVu3buH8+fMwNTWt8piioqKwfPly3L17F23atCmzbm5uLgYMGICPPvoIL730EtavX493330XFhYWSEtL0/mYixcvRnW4pdw///wDLy8vfPPNNxg1ahTWrFmD/v37a74HCwoK9NqfcSXFSWX44IMP0LJlS/zyyy9QKBRa6xITEyWKqnwFBQVQq9WSfOh1YWxsDGPjqnlL7969GytXrsSwYcOwc+dOyOVyzTp/f3/06dMHr7/+On7//XetmIyNjeHr66tZnjJlCp577jns3r0bQUFBcHJyqpL4iYDC5FhXT35XVbX+/fvj8OHDePDgAezt7TXlu3btgpOTE5o0aYKUlJQqjalTp054+PAh7OzscODAAbz++uul1l21ahV++uknnD17FiqVqkLHa9++PSIiIvDNN99gyJAhFQ3bIA4fPozY2FgcOXIEAwYM0JTb2dlh8eLFuHz5Mjp06KDz/mply8jdu3cxZswYODk5QalUolWrVti8ebNWHbVajYULF6JTp06wtraGhYUFevbsiR9//LHY/vbs2YNOnTrB0tISVlZWaNOmDT799FMAhc3URW/APn36aJrSymqq+/vvv9GlS5cSP9yOjo5ay6mpqfD394e1tTVsbGzg5+eHiIgIyGQybN26VVOvtP7ckvrUV6xYgeeeew5169aFmZkZOnXqhAMHDhTbViaTYdq0adi5cydatWoFpVKJEydOANDtHAPAmjVr0KpVK5ibm8PW1hadO3cusytKCAF7e3sEBARoygoKCmBjYwO5XI7U1FRN+fLly2FsbIx//vkHQPExIzKZDJmZmdi2bZvmdXmyf7zo/NrY2MDa2hqjR49GVlZWqfEVWbRoEWxtbfHll19qJSIAoFKpMGfOHFy+fBmHDh0qcz8ymQw9evSAEALR0dHlHve7777DgAED4OrqCqVSiUaNGuHDDz9Efn6+Vr3evXujdevWuHbtGvr06QNzc3PUq1cPH3/8cbF93rlzB6+++iosLCzg6OiIWbNmIScnp9xYqtq3336L1q1bw9TUFK1bt8Y333xTYr2CggKsXr0arVq10jSFT5w4scQ/VMePH0evXr00n+0uXboUe3/u378fnTp1gpmZGezt7eHr64u7d+9q1fH390edOnUQFxeHV155BXXq1EG9evWwbt06AMAff/yBF154ARYWFmjYsGGxYxR1d/3888+YOHEi6tatCysrK4waNarEuNevX6/5TLq6umLq1Klanw0AuHnzJl577TU4OzvD1NQU9evXxxtvvKH1K/zfY0bK+y4r6TsmMTERY8eOhZOTE0xNTdGuXTts27ZNq86/u+q+/PJLNGrUCEqlEl26dMFvv/1W7LmVZvDgwVAqldi/f79W+a5duzB06NBin8OqYGlpCTs7u3LrFRQU4NNPP8V//vMfqFQq5OXl6fQ986Q33ngDTZs2rRatI+np6QBQ7AeUi4sLAMDMzEyv/dW6lpGEhAR069ZN84fUwcEBx48fx9ixY5Geno633noLQOGJ/Oqrr/Dmm29i/PjxyMjIwKZNm+Dj44OwsDC0b98eAHD69Gm8+eab6Nu3L5YvXw4AuH79Os6dO4eZM2fi+eefx4wZM/DZZ5/hvffeQ4sWLQBA829JGjZsiJCQENy5cwf169cvtZ4QAoMHD8bZs2cxadIktGjRAt988w38/Pye6hx9+umnGDRoEEaMGAG1Wo09e/bg9ddfL5bhAsAPP/yAffv2Ydq0abC3t4e7u7vO53jjxo2YMWMG/vvf/2LmzJl49OgRrly5gl9//RXDhw8vMTaZTIbu3bvj559/1pRduXIFaWlpMDIywrlz5zQxhoaGokOHDqhTp06J+9q+fTvGjRsHlUqFCRMmAAAaNWqkVWfo0KHw8PDA0qVLER4ejq+++gqOjo6a17okN2/eRFRUFPz9/WFlZVVinVGjRiEwMBDff/89hg4dWuq+AGiaxG1tbcusBxT+wahTpw4CAgJQp04d/PDDD1i4cCHS09PxySefaNVNSUnBSy+9hCFDhmDo0KE4cOAA5syZgzZt2uDll18GAGRnZ6Nv376Ii4vDjBkz4Orqiu3bt+OHH34oN5aqdOrUKbz22mto2bIlli5diocPH2L06NElfn4mTpyIrVu3YvTo0ZgxYwZiYmKwdu1aXLp0CefOnYOJiQmAwnM5ZswYtGrVCvPmzYONjQ0uXbqEEydOaN6fRfvp0qULli5dioSEBHz66ac4d+4cLl26BBsbG81x8/Pz8fLLL+P555/Hxx9/jJ07d2LatGmwsLDA/PnzMWLECAwZMgQbNmzAqFGj4OXlBQ8PD63Yp02bBhsbG3zwwQeIiorC559/jtjYWM24CaAw6V60aBG8vb0xefJkTb3ffvtN8/zUajV8fHyQk5OD6dOnw9nZGXfv3sWRI0eQmpoKa2vrYudN3++y7Oxs9O7dG3/99RemTZsGDw8P7N+/H/7+/khNTcXMmTO16u/atQsZGRmYOHEiZDIZPv74YwwZMgTR0dGa16Qs5ubmGDx4MHbv3o3JkycDAC5fvow///wTX331Fa5cuVLuPgAgKytLp0RALpfr9JnUxbVr13Dv3j20bdsWEyZMwLZt26BWqzU/bPv06aPTfuRyORYsWIBRo0ZVqHUkNzdX5y4hOzs7GBmV3l7x/PPPw8jICDNnzsTKlStRv359XLlyBR999BFeffVVNG/eXK/YIGqQLVu2CADit99+K7XO2LFjhYuLi3jw4IFW+RtvvCGsra1FVlaWEEKIvLw8kZOTo1UnJSVFODk5iTFjxmjKZs6cKaysrEReXl6px9y/f78AIH788UednsemTZsEAKFQKESfPn3E+++/L0JDQ0V+fr5WvW+//VYAEB9//LGmLC8vT/Ts2VMAEFu2bNGU9+rVS/Tq1avYsfz8/ETDhg21yorOQRG1Wi1at24tXnjhBa1yAMLIyEj8+eefWuW6nuPBgweLVq1alXkuSvLJJ58IuVwu0tPThRBCfPbZZ6Jhw4ZCpVKJOXPmCCGEyM/PFzY2NmLWrFma7QIDA8WTb2kLCwvh5+dX7BhFdf/9WgshxH/+8x9Rt27dMuMrel1WrVpVZj0rKyvRsWNHzbKfn5+wsLAQSUlJIikpSfz1119ixYoVQiaTidatW4uCgoIy9ydE8ddOCCEmTpwozM3NxaNHjzRlvXr1EgDE119/rSnLyckRzs7O4rXXXtOUrV69WgAQ+/bt05RlZmaKxo0b6/SeLu8zGRMTIwCITz75pMT1n3zyiQAgYmJiyjxO+/bthYuLi0hNTdWUnTp1SgDQen+HhoYKAGLnzp1a2584cUKrPDU1VVhaWoquXbuK7OxsrbpFr4NarRaOjo6idevWWnWOHDkiAIiFCxdqyvz8/AQAsWTJEk1ZSkqKMDMzEzKZTOzZs0dTHhkZKQCIwMBATVnReezUqZNQq9Wa8o8//lgAEN99950QQojExEShUChEv379tL4v1q5dKwCIzZs3CyGEuHTpkgAg9u/fX+Z5bdiwodbno6zvsie/Y4reOzt27NCUqdVq4eXlJerUqaP5/Ba9B+rWrSuSk5M1db/77jsBQHz//fdlxvjjjz9qnsuRI0eETCYTcXFxQggh3nnnHeHp6amJT5fvm6LPfnmPJ783y1PWuTt06JDmHDRp0kRs2bJFbNmyRTRp0kQoFApx+fLlMvf9789RXl6eaNKkiWjXrp3mvVr0nJKSksrcT9G51OVR3mdSCCG++uorYWNjo7Wdn5+fyM3NLXfbJ9WqbhohBA4ePIiBAwdCCIEHDx5oHj4+PkhLS0N4eDiAwgyzqJukoKAAycnJyMvLQ+fOnTV1AMDGxgaZmZk4ffq0weIcM2YMTpw4gd69e+Ps2bP48MMP0bNnTzRp0gTnz5/X1Dt27BiMjY01vwKK4p4+ffpTHf/fzWcpKSlIS0tDz549tZ53kV69eqFly5aaZX3OsY2NDe7cuaNXUywA9OzZE/n5+ZpzERoaip49e6Jnz54IDQ0FAFy9ehWpqano2bOn3s//3yZNmlTs2A8fPtQ0QZYkIyMDQGETbVksLS01dYtkZmbCwcEBDg4OaNy4MWbPno3u3bvju+++0+my5H+/dhkZGXjw4AF69uyJrKwsREZGatWtU6eO1vgUhUIBlUql1R107NgxuLi44L///a+mzNzcXNOSVB3cv38fERER8PPz0/pF/+KLL2q9N4HCLhVra2u8+OKLWu/NTp06oU6dOppu2NOnTyMjIwNz584tNgaq6HX4/fffkZiYiClTpmjVGTBgAJo3b46jR48Wi3XcuHGa/9vY2KBZs2awsLDQah1r1qwZbGxsSuyWmzBhglYrweTJk2FsbIxjx44BAIKDg6FWq/HWW29p/WodP348rKysNDEVnaeTJ09WqDtAF8eOHYOzszPefPNNTZmJiQlmzJiBf/75Bz/99JNW/WHDhmm1NBR9dnXpnizSr18/2NnZYc+ePRBCYM+ePVrH18WoUaNw+vTpch87d+7Ua79lKepKzsjIQEhICPz9/eHv74/g4GAIIUrsPi1NUevI5cuX9R6c365dO52e++nTp+Hs7Fzu/urVqweVSoXVq1fjm2++QUBAAHbu3Im5c+fqFRdQy7ppkpKSkJqaii+//BJffvlliXX+PUB027ZtWLlyJSIjI5Gbm6sp/3fT6ZQpU7Bv3z7NZbj9+vXD0KFD8dJLLz1VrD4+PvDx8UFWVhYuXryIvXv3YsOGDXjllVcQGRkJR0dHxMbGwsXFpVg3RLNmzZ7q2EeOHMH//vc/REREaI0NKOmP4ZPNyPqc4zlz5iA4OBgqlQqNGzdGv379MHz4cHTv3r3M+Dp27Ahzc3OEhobCx8cHoaGhWLRoEZydnbFmzRo8evRIk5T06NFDr+f+JDc3N63loi/LlJSUUrtgipKQJxONJ2VkZBQbr2Nqaorvv/8eQOFYjY8//hiJiYk696/++eefWLBgAX744YdiCdOTza/169cv9pra2tpqNWfHxsaicePGxeo97XtMX2UlYrGxsQCAJk2aFFvXrFkzrST65s2bSEtLKzb2qkjRe/Pvv/8GgDIvBS06bknnonnz5jh79qxWmampKRwcHLTKrK2tS3wdrK2tSxwL8uRzrFOnDlxcXDRdeaXFpFAo4OnpqVnv4eGBgIAABAUFYefOnejZsycGDRoEX1/fErtoKiI2NhZNmjQp1pRf1K1TFEuRsj5rujIxMcHrr7+OXbt2QaVS4fbt26V2+ZbG09MTnp6eem3ztIo+3927d0eDBg005W5ubujRo4fWj1BdjBgxAh9++CEWL16MV199VeftbG1t4e3trdexSnPu3Dm88sor+OWXX9C5c2cAwKuvvgorKyssWrQIY8aMKfZjoSy1KhkpupTI19e31HEVbdu2BQDs2LED/v7+ePXVV/HOO+/A0dERcrkcS5cu1XxRAYUDSiMiInDy5EkcP34cx48fx5YtWzBq1KhiA7UqwtzcXPOr397eHosWLcLx48f1Hhcik8lKHND05MDG0NBQDBo0CM8//zzWr18PFxcXmJiYYMuWLSUOLH3yj6Q+57hFixaIiorCkSNHcOLECRw8eBDr16/HwoULsWjRolKfi4mJCbp27Yqff/4Zf/31F+Lj49GzZ084OTkhNzcXv/76K0JDQ9G8efNiX/76Km3QW0nnskjRB6ysPurY2Fikp6cX+9KTy+VaXwY+Pj5o3rw5Jk6ciMOHD5cZa2pqKnr16gUrKyssXrwYjRo1gqmpKcLDwzFnzpxil9JV5LkZWlGLwpOXXBcp+tVuqCu0CgoK4OjoWOqv2qd9v5SltPMt1euwcuVK+Pv747vvvsOpU6cwY8YMLF26FL/88kuZY9Uqi6HOw/Dhw7FhwwZ88MEHaNeunV5/8IDCVoqiloqyyOVyg71fXF1dARQf7AkU/o25dOmSXvsrah0pen11pVarS5wDqSQODg5lDgr+4osv4OTkpElEigwaNAgffPABzp8//+wmIw4ODrC0tER+fn652d+BAwfg6emJQ4cOaf1qCQwMLFZXoVBg4MCBGDhwIAoKCjBlyhR88cUXeP/990v8VVlRRS/q/fv3Afz/QNd//vlHq3UkKiqq2La2trYlNnc++evk4MGDMDU1xcmTJ7Uu69uyZYtOMepzjgHAwsICw4YNw7Bhw6BWqzFkyBB89NFHmDdvXpl/gHr27Inly5cjODgY9vb2aN68OWQyGVq1aoXQ0FCEhobilVdeKff4hnpt/q1JkyZo1qwZvv32W3z66acldtd8/fXXAFDmpX5A4cjzWbNmYdGiRfjll1/QrVu3UuueOXMGDx8+xKFDh/D8889rymNiYir4TArfY1evXoUQQutclfQeqwgHBweYm5uXur+oqCiYm5trXapZUowASpy35cn9NmrUCMHBwejevXuZrU1FA5mvXr2Kxo0bl3ncqKgovPDCC8WOW7TekG7evKk1mPGff/7B/fv30b9//2Ix/TvRVavViImJKfaZbNOmDdq0aYMFCxbg/Pnz6N69OzZs2ID//e9/JR5fn89Lw4YNceXKFRQUFGi1jhR1F1bG+QEKW0Pd3Nxw5syZMgeal2bFihVl/hgq0rBhQ4NN/tamTRuYmJgUuwoLAO7du1ehpMfX1xf/+9//sGjRIgwaNEinbc6fP6/zYNmYmJgyZzdOSEgo9mMXgKaXIS8vT6fjFKlVY0bkcjlee+01HDx4EFevXi22PikpSasuoJ2V//rrr7hw4YLWNk/OJGdkZKT55V/UxWFhYQEAxS6tK01ISEiJ5UX9wkVNsP3790deXh4+//xzTZ38/HysWbOm2LaNGjVCZGSk1nO8fPkyzp07p1VPLpdDJpNpvYlu3bqlc9+jPuf4yXOnUCjQsmVLCCG0usVK0rNnT+Tk5GD16tXo0aOH5kuyZ8+e2L59O+7du6fTeBELCwudXxd9BAYGIiUlBZMmTSr2gbx48SKWL1+ODh06aK5aKcv06dNhbm6OZcuWlVmvpPesWq3G+vXrK/AMCvXv3x/37t3TurQ7Kyur1C44fcnlcvTr1w/ff/894uLitNbFxcXh+++/R79+/cr8Bebi4oL27dtj27ZtWl1Rp0+fxrVr17TqDh06FPn5+fjwww+L7ScvL0/zXujXrx8sLS2xdOnSYrPfFp3fzp07w9HRERs2bNDqzjx+/DiuX79e7MozQ/jyyy+1Phuff/458vLyNO8jb29vKBQKfPbZZ1rvg02bNiEtLU0TU3p6erE/Bm3atIGRkVGZl23r813Wv39/xMfHY+/evZqyvLw8rFmzBnXq1EGvXr3Kf8IVIJPJ8NlnnyEwMBAjR47Ue3spxoxYWlqif//+OH/+vNbYruvXr+P8+fN48cUX9d5nUetIREREua2qRQw5ZqRp06ZISEgoNo3F7t27AUCvOUaAGtoysnnzZs18F/82c+ZMLFu2DD/++CO6du2K8ePHo2XLlkhOTkZ4eDiCg4M1TVSvvPIKDh06hP/85z8YMGAAYmJisGHDBrRs2VKrCW/cuHFITk7GCy+8gPr16yM2NlYzw2ZR32j79u0hl8uxfPlypKWlQalU4oUXXii133rw4MHw8PDAwIED0ahRI2RmZiI4OBjff/89unTpgoEDBwIABg4ciO7du2Pu3Lm4desWWrZsiUOHDpV4adaYMWMQFBQEHx8fjB07FomJidiwYQNatWqlNbZgwIABCAoKwksvvYThw4cjMTER69atQ+PGjXW+NE7Xc9yvXz84Ozuje/fucHJywvXr17F27VoMGDCg3MGfXl5eMDY2RlRUlNZgyueff16TnOmSjHTq1AnBwcEICgqCq6srPDw80LVrV52eZ1nefPNN/P777wgKCsK1a9cwYsQI2NraIjw8HJs3b4aDgwMOHDig0yRsdevWxejRo7F+/Xpcv3691Espn3vuOdja2sLPzw8zZsyATCbD9u3bn6q5f/z48Vi7di1GjRqFixcvwsXFBdu3b4e5uble+ynrM7lkyRJ069YNHTt2xIQJE+Du7o5bt27hyy+/hEwmw5IlS8rd/9KlSzFgwAD06NEDY8aMQXJysmYOm39/Xnv16oWJEydi6dKliIiIQL9+/WBiYoKbN29i//79+PTTT/Hf//4XVlZWWLVqFcaNG4cuXbpg+PDhsLW1xeXLl5GVlYVt27bBxMQEy5cvx+jRo9GrVy+8+eabmkt73d3dMWvWLL3OkS7UajX69u2LoUOHIioqCuvXr0ePHj00v3wdHBwwb948LFq0CC+99BIGDRqkqdelSxfNgOUffvgB06ZNw+uvv46mTZsiLy8P27dv1/yYKI0+32UTJkzAF198AX9/f1y8eBHu7u44cOAAzp07h9WrV5f7GX8agwcPxuDBgyu0raHHjBS1Mv35558ACqcUKBpPtGDBAk29JUuWICQkBC+88AJmzJgBAPjss89gZ2eH9957r0LHLho7EhERoVN9Q44ZmTZtGrZs2YKBAwdi+vTpaNiwIX766Sfs3r0bL774ov7fs3pffyOhosvfSnvcvn1bCCFEQkKCmDp1qmjQoIEwMTERzs7Oom/fvuLLL7/U7KugoEAsWbJENGzYUCiVStGhQwdx5MiRYpfCHjhwQPTr1084OjoKhUIh3NzcxMSJE8X9+/e1Ytu4caPw9PQUcrm83Esid+/eLd544w3RqFEjYWZmJkxNTUXLli3F/PnzNZfDFXn48KEYOXKksLKyEtbW1mLkyJGay/b+fWmvEELs2LFDeHp6CoVCIdq3by9OnjxZ4qW9mzZtEk2aNBFKpVI0b95cbNmypcTLYgGIqVOnlvgcdDnHX3zxhXj++edF3bp1hVKpFI0aNRLvvPOOSEtLK/Xc/FuXLl0EAPHrr79qyu7cuSMAiAYNGhSrX9JziIyMFM8//7wwMzPTXHb277pPXgpX9B7T5bI2IYQ4fPiw8Pb21rq8rVWrViU+x6JLe0vy999/C7lcXuJlyP927tw50a1bN2FmZiZcXV3Fu+++K06ePFnsPVfaZY4lvR9iY2PFoEGDhLm5ubC3txczZ87UXAqr66W95X0mr1+/LoYNGyYcHR2FsbGxcHR0FG+88Ya4fv16mfv/t4MHD4oWLVoIpVIpWrZsKQ4dOlTi8xFCiC+//FJ06tRJmJmZCUtLS9GmTRvx7rvvinv37mnVO3z4sHjuueeEmZmZsLKyEiqVSuzevVurzt69e0WHDh2EUqkUdnZ2YsSIEeLOnTtadUp7bUt7HRo2bCgGDBhQ7Dz+9NNPYsKECcLW1lbUqVNHjBgxQjx8+LDY9mvXrhXNmzcXJiYmwsnJSUyePFmkpKRo1kdHR4sxY8aIRo0aCVNTU2FnZyf69OkjgoODi8Xx5HuutO+ykqYPSEhIEKNHjxb29vZCoVCINm3aFPteKuvybjxxiXNJ/n1pb1l0vbTX0Mp6/z/p4sWLwtvbW1hYWAhLS0sxePBgcePGjXKPUdY5/PdnsLxLew0tMjJS/Pe//9X8HWjYsKGYPXu2yMzM1HtfMiGqwST3pJdbt27Bw8MDW7Zs0fmOm1Q1xo0bh02bNmHjxo1al3kSlaVocrXffvut2IBAomdBjeymIaquvvjiCyQkJGDy5MlwdXXVDDwkIqLSMRkhMiC5XK6ZR4SIiHRTq66mISIiopqHY0aIiIhIUmwZISIiIkkxGSEiIiJJ1YgBrAUFBbh37x4sLS0rZXpvIiIiMjwhBDIyMuDq6lrspor/ViOSkXv37mnd6ZCIiIhqjtu3b5d5g8YakYwUTSt8+/btUm/rTkRERNVLeno6GjRoUO7tAWpEMlLUNWNlZcVkhIiIqIYpb4gFB7ASERGRpJiMEBERkaSYjBAREZGkasSYESIiqp2EEMjLy0N+fr7UoVAFyOVyGBsbP/W0G0xGiIhIEmq1Gvfv30dWVpbUodBTMDc3h4uLCxQKRYX3wWSEiIiqXEFBAWJiYiCXy+Hq6gqFQsFJLWsYIQTUajWSkpIQExODJk2alDmxWVmYjBARUZVTq9UoKChAgwYNYG5uLnU4VEFmZmYwMTFBbGws1Go1TE1NK7QfDmAlIiLJVPSXNFUfhngN2TJSxfIL8hGeGI6krCQ4mDugo2NHyI3kUodFREQkmWc2GVFnP0LYhh3ISkyFuaMNVJN8oTCrWPOSroJjg7EsbBkSshI0ZU7mTpirmgvvht6VemwiIqLqSu+2lZ9//hkDBw6Eq6srZDIZvv3223K3OXPmDDp27AilUonGjRtj69atFQjVcEICV+HrKUdxOcoTN1M64nKUJ76echQhgasq7ZjBscEIOBOglYgAQGJWIgLOBCA4NrjSjk1EVFvlFwhc+Pshvou4iwt/P0R+gaj0YwohMGHCBNjZ2UEmkyEiIgIPHz6Eo6Mjbt26pdM+1Go13N3d8fvvv1dusDWE3slIZmYm2rVrh3Xr1ulUPyYmBgMGDECfPn0QERGBt956C+PGjcPJkyf1DtYQQgJXITK+LXIUNlrlOQobRMa3rZSEJL8gH8vClkFAwCgf8Lrhif4RHeB1wxOy/MIPzvKw5cgv4HX2RES6OnH1Pnos/wFvbvwFM/dE4M2Nv6DH8h9w4ur9yj3uiRPYunUrjhw5gvv376N169b46KOPMHjwYLi7u+u0D4VCgdmzZ2POnDmVGmtNoXc3zcsvv4yXX35Z5/obNmyAh4cHVq5cCQBo0aIFzp49i1WrVsHHx0ffwz8VdfYjRMe5AQoAT15CJpMBQiA6rgF6Zj8yaJdNeGI4ErIS0O+PNmiW/BpyFbYAALdsoPP9FETZHcSpNn8gPDEcXZy7GOy4RES11Ymr9zF5RziebAeJT3uEyTvC8blvR7zU2qVSjv3333/DxcUFzz33HAAgKysLmzZt0vtH9ogRI/D222/jzz//RKtWrSoj1Bqj0ocxX7hwAd7e2uMhfHx8cOHChVK3ycnJQXp6utbDEMI+3w610rZ4IlJEJoNaaYewz7cb5HhFkrKS0O+PNvDMGItcExutdbkmNvDMGIt+f7RBUlaSQY9LRFQb5RcILPr+WrFEBICmbNH31yqly8bf3x/Tp09HXFwcZDIZ3N3dcezYMSiVSnTr1k1Tb/HixXB1dcXDhw81ZUW9BAUFBQAAW1tbdO/eHXv27DF4nDVNpScj8fHxcHJy0ipzcnJCeno6srOzS9xm6dKlsLa21jwaNGhgkFjuxt81aD1d2RpZo1nya4ULJbXIAGiWPAS2RtYGPS4RUW0UFpOM+2mPSl0vANxPe4SwmGSDH/vTTz/F4sWLUb9+fdy/fx+//fYbQkND0alTJ6168+fPh7u7O8aNGwcAWLduHc6fP49t27ZpXQqrUqkQGhpq8Dhrmmp5gfe8efOQlpamedy+fdsg+1Xb6Ta7n671dJV/6Fph10wZLTK5CjvkH7pm0OMSEdVGiRmlJyIVqacPa2trWFpaQi6Xw9nZGQ4ODoiNjYWrq6tWPblcjh07diAkJARz587FO++8g3Xr1sHNzU2rnqurK2JjYw0eZ01T6cmIs7MzEhK0ryBJSEiAlZUVzMzMStxGqVTCyspK62EIjd/whok6BRClNN0JARN1Mhq/YdjLbLOS0gxaj4joWeZoqduYPl3rPa3s7OwSZx719PTEihUrsHz5cgwaNAjDhw8vVsfMzIz35kEVJCNeXl4ICQnRKjt9+jS8vLwq+9DFqNy6IdbucOHCkwnJ4+VYuyNQuXWDIZk72hi0HhHRs0zlYQcXa1OU1oYtA+BibQqVh12VxGNvb4+UlJQS1/3888+Qy+W4desW8vLyiq1PTk6Gg4NDZYdY7emdjPzzzz+IiIhAREQEgMJLdyMiIhAXFwegsItl1KhRmvqTJk1CdHQ03n33XURGRmL9+vXYt28fZs2aZZhnoAe5kRx9po9EtOUmmOSmaq0zyU1BtOUm9Jnua/AZUVWTfKHMKbtFRpmTDNUkX4Mel4ioNpIbyRA4sCUAFEtIipYDB7aE3KhqbrzXoUMHXLtWvJt97969OHToEM6cOYO4uDh8+OGHxepcvXoVHTp0qIowqzW9k5Hff/8dHTp00Jy8gIAAdOjQAQsXLgQA3L9/X5OYAICHhweOHj2K06dPo127dli5ciW++uqrKr+st4h3Q2+8PH0M9vf7HJfrfoo4s624XPdT7O+3AS9PH1MpM6EqzEzh4fb4nJTSIuPhdrvSZ4AlIqotXmrtgs99O8LZWvt709natFIv6y2Jj48P/vzzT63WkTt37mDy5MlYvnw5evTogS1btmDJkiX45ZdftLYNDQ1Fv379qizW6komRGk/16uP9PR0WFtbIy0tzWDjR6S4R0xI4CrExLkhR2mrKVPmJMPD7Tb6LqrcliIppr8nIirNo0ePEBMTAw8Pjwrf6RUovMw3LCYZiRmP4GhZ2DVT2S0iq1evxurVq7VmW+3atSvGjBmDiRMnQgiBF198EcbGxjh+/Dhkjy9emDFjBo4dO4aIiAjUqVMHFy5cQP/+/XHv3r1Sx1DWBGW9lrr+/X5mkxGpSJEUlJwEpcDDLa7SkyAiopIYKhmpLo4ePYp33nkHV69e1fkutsOGDUO7du3w3nvvVXJ0lcsQycgze6M8qSjMTNFj1rgqO17R9PdQaJcXTn9vAwSuYkJCRPSUBgwYgJs3b+Lu3bs6zY2lVqvRpk0bScZPVkdMRmoxqaa/JyJ6Fr311ls611UoFFiwYEHlBVPDVMtJz8gwpJr+noiISB9MRmoxqaa/JyIi0geTkVpMqunviYiI9MFkpBaTavp7IiIifTAZqcWkmv6eiIhIH0xGajGppr8nIiLSBy/treW8G3oD04HlF5bB4w8L2GZZI8U8DTFtMjHHa26lTH9PRESkDyYjzwDvht7o06APwntX7fT3RESVriAfiD0P/JMA1HECGj4H1KDvtpEjR6JFixY6z8I6d+5cZGZmYs2aNZUcWdViN80zQm4kRxfnLujv2R9dnLswESGimu/aYWB1a2DbK8DBsYX/rm5dWC6RNm3aYNKkSSWu2759O5RKJR48eAAAuHz5Mo4dO4YZM2bovP/Zs2dj27ZtiI6O1iuubdu2oUuXLjA3N4elpSV69eqFI0eOaNU5c+YMZDKZ5uHg4ID+/fvjjz/+0OtYFcFkhCqVOvsRzq76CqfmrcDZVV9Bnf1I6pCIqDa4dhjYNwpIv6ddnn6/sFyihGTs2LHYs2cPsrOzi63bsmULBg0aBHt7ewDAmjVr8Prrr6NOnTo679/e3h4+Pj74/PPPdd5m9uzZmDhxIoYNG4YrV64gLCwMPXr0wODBg7F27dpi9aOionD//n2cPHkSOTk5GDBgANRqtc7HqwgmI1RpQgJX4espR3E5yhM3UzricpQnvp5yFCGBq6QOjYhqsoJ84MQcACVNW/C47MTcwnqVICMjAyNGjICFhQVcXFywatUq9O7dG2+99RZ8fX2RnZ2NgwcPam0TExODM2fOYOzYsQCA/Px8HDhwAAMHDtTUiYyMhLm5OXbt2qUp27dvH8zMzHDt2jVN2cCBA7Fnzx6dYv3ll1+wcuVKfPLJJ5g9ezYaN26MFi1a4KOPPsJbb72FgIAA3L59W2sbR0dHODs7o2PHjnjrrbdw+/ZtREZG6n2e9MFkhCpF0Q36chQ2WuWFN+hry4SEiCou9nzxFhEtAki/W1ivEgQEBODcuXM4fPgwTp8+jdDQUISHhwMobLkYPHgwNm/erLXN1q1bUb9+ffTr1w8AcOXKFaSlpaFz586aOs2bN8eKFSswZcoUxMXF4c6dO5g0aRKWL1+Oli1bauqpVCrcuXMHt27dKjfW3bt3o06dOpg4cWKxdW+//TZyc3OLJU5F0tLSNEmPQqEosY6hcAArGZw6+xFiyrlBX0xcA6h5gz4iqoh/EgxbTw8ZGRnYtm0bdu3ahb59+wIo7H5xdXXV1Bk7dixefvllxMTEwMPDA0IIbNu2DX5+fjAyKmwDiI2NhVwuh6Ojo9b+p0yZgmPHjsHX1xcKhQJdunTB9OnTteoUHSs2Nhbu7u5lxnvjxg00atSoxGTC1dUVVlZWuHHjhlZ5/fr1AQCZmZkAgEGDBqF58+blnZqnwpYRMriwDTuQU84N+nKUdgjbsKNqAyOi2qGOk2Hr6SE6Ohq5ublQqVSaMmtrazRr1kyz/OKLL6J+/frYsmULACAkJARxcXEYPXq0pk52djaUSiVkJXxPbt68GVeuXEF4eDi2bt1arI6ZmRkAICsrS6eYRWmzcD/2ZKISGhqKixcvYuvWrWjatCk2bNig03GeBpMRMrishBSD1iMi0tLwOcDKFUBp99WSAVb1CutJwMjICP7+/ti2bRsKCgqwZcsW9OnTB56enpo69vb2yMrKKnFg6OXLl5GZmYnMzEzcv3+/2Prk5GQAgIODQ7mxNGnSBNHR0SUe5969e0hPT0fTpk21yj08PNCsWTP4+flh3LhxGDZsWLnHeVpMRsjg0u10GzSmaz0iIi1GcuCl5Y8XnkxIHi+/tKxS5hvx9PSEiYkJfvvtN01ZWlpasa6O0aNH4/bt2zh06BC++eYbzcDVIu3btwcArYGpQGGi4e/vj/nz58Pf3x8jRowodmXO1atXYWJiglatWpUb75tvvol//vkHX3zxRbF1K1asgKmpaZnJxtSpU3H16lV888035R7raTAZIYNTDmyh0w36lANbVG1gRFR7tBwEDP0asHLRLrdyLSxvOahSDmtpaQk/Pz+88847+PHHH/Hnn39i7NixMDIy0upO8fDwwAsvvIAJEyZAqVRiyJAhWvtxcHBAx44dcfbsWa3ySZMmoUGDBliwYAGCgoKQn5+P2bNna9UJDQ1Fz549Nd01ZfHy8sLMmTPxzjvvYOXKlfj7778RGRmJBQsW4LPPPsPGjRtRt27dUrc3NzfH+PHjERgYWG53z9NgMkIG52zriii7x6OzS7lBX5TdITjbuoKIqMJaDgLeugr4HQFe21T471t/VFoiUiQoKAheXl545ZVX4O3tje7du6NFixYwNdUekD927FikpKRg+PDhxdYBwLhx47Bz507N8tdff41jx45h+/btMDY2hoWFBXbs2IGNGzfi+PHjmnp79uzB+PHjdY539erVWL9+PXbv3o3WrVujRYsW+OSTT/DDDz/A19e33O2nTZuG69evY//+/TofU18yUZmpjoGkp6fD2toaaWlpsLKykjocKkd+QT58Dvqg3a+OaJb8GnIVtpp1JupkRNkdwpWuSTjx2gnOBEv0jHr06JHmapOS/lDXJJmZmahXrx5WrlxZrDumLNnZ2WjWrBn27t0LLy8vnbY5fvw43n77bVy5cgXGxhW7IPbWrVvo1asXvLy8sHPnTsjlT/c9XNZrqevfb17aSwYnN5JjrmouArICEJJ/Faq/PTQ36AtrFIMCORCkCmIiQkQ10qVLlxAZGQmVSoW0tDQsXrwYADB48GC99mNmZoavv/5aMz28LjIzM7Fly5YKJyIA4O7ujjNnzmDbtm2IiIhAp06dKrwvQ2HLCFWa4NhgLAtbhoSs/7/W39ncGXNUc3i3YKJnXE1uGbl06RLGjRuHqKgoKBQKdOrUCUFBQWjTpk2VxzJp0iTs2FHyNAm+vr5VclmuIVpGmIxQpcovyEd4Iu8WTETaanIyUp0kJiYiPT29xHVWVlbFJlWrDOymoWqv6G7BRERkeI6OjlWScFQ2Xk1DREREkmIyQkRERJJiMkJERESS4pgRqrU4eJaIqGZgMkK1UnBsMJZfWAaPPyw0c5zEtMnEHK+5vKyYiKiaYTJCtU5wbDCOr9mM15Mna2Z/dcsGOp9KwfHfNwPTwYSEqJaoqS2gUVFR6NWrF27evAlLS8ty6z948AAtW7ZEeHg46tevXwURVi2OGaFaJb8gHz+u2Q7PjLHINbHRWpdrYgPPjLH4cc0O5BfwjsFENV1wbDB8DvpgzMkxmBM6B2NOjoHPQR8ExwZLGlfv3r0hk8mKPfLy8jR15s2bh+nTp+uUiACAvb09Ro0ahcDAQL1iuX37NsaMGQNXV1coFAo0bNgQM2fOxMOHD0uN2dTUFE2bNsXSpUsr9eZ4/8ZkhGqVsLhf0DD58U2yZE/cWvzxcsPkVxAW90sVR0ZEhhQcG4yAMwFaMzwDQGJWIgLOBEiekIwfPx7379/XehRN4R4XF4cjR47A399fr32OHj0aO3fuRHJysk71o6Oj0blzZ9y8eRO7d+/GX3/9hQ0bNiAkJAReXl7F9lMUc1RUFObNm4eFCxdWyQyuAJMRqmX+2hNc2DXzZCJSRCZDrsIOf+2R9ouKiCouvyAfy8KWQaD4r/aisuVhyyutBTQjIwMjRoyAhYUFXFxcsGrVKvTu3RtvvfWWpo65uTmcnZ21HkX27duHdu3aoV69epqyMWPGoG3btsjJyQEAqNVqdOjQAaNGjdLUadWqFVxdXfHNN9/oFOfUqVOhUChw6tQp9OrVC25ubnj55ZcRHByMu3fvYv78+Vr1i2Ju2LAhRo8ejbZt2+L06dMVOUV6YzJCtYoiWbcmRV3rEVH1E54YXqxF5N8EBOKz4hGeGF4pxw8ICMC5c+dw+PBhnD59GqGhoQgP1/1YoaGh6Ny5s1bZZ599hszMTMydOxcAMH/+fKSmpmLt2rVa9VQqFUJDQ8s9RnJyMk6ePIkpU6bAzMxMa52zszNGjBiBvXv3ltgNI4RAaGgoIiMjoVAodH5eT4PJCNUq9ZzrlV9Jj3pEVP0kZSUZtJ4+MjIysG3bNqxYsQJ9+/ZF69atsWXLFuTna7fCrF+/HnXq1NE83n77bc262NhYuLq6atWvU6cOduzYgXXr1mHhwoVYvXo1tm/fXux+Lq6uroiNjS03zps3b0IIgRYtWpS4vkWLFkhJSUFS0v+fo6KYlUolnn/+eRQUFGDGjBnlHssQeDUN1SqqySNxfcpRqBU2JXfVCAGFOgWqySOrPDYiMgwHcweD1tNHdHQ0cnNzoVKpNGXW1tZo1qyZVr0RI0ZodYPY2Nho/p+dnV3izQG9vLwwe/ZsfPjhh5gzZw569OhRrI6ZmRmysrJ0jre8Aaj/bvkoijklJQWBgYF47rnn8Nxzz+l8rKfBZIRqFYWZKTzd4hAZbwMIoZ2QPP5QerrdhsKMdwklqqk6OnaEk7kTErMSSxw3IoMMTuZO6OjYUYLoCllbW6Nx48YlrrO3t0dKSkqx8oKCApw7dw5yuRx//fVXidsmJyfDwaH8JKtx48aQyWS4fv06/vOf/xRbf/36dTg4OGglSf+Oed++fWjcuDG6desGb+/KnwqB3TRU6/RdNAvNna9AqU7VKleqU9Dc+Qr6LpolTWBEZBByIznmqgrHVsig3QJatDxHNadS5hvx9PSEiYkJfvvtN01ZWloabty4ofM+OnTogGvXrhUr/+STTxAZGYmffvoJJ06cwJYtW4rVuXr1Kjp06FDuMerWrYsXX3wR69evR3Z2tta6+Ph47Ny5s8yreerUqYOZM2di9uzZVXJ5L5MRqpX6LpqFUesHoF2zaDSxDUe7ZtEYtf4VJiJEtYR3Q28E9Q6Co7mjVrmTuROCegdV2sSGlpaW8PPzwzvvvIMff/wRf/75J8aOHQsjIyPISruK7wk+Pj64cOGC1jiTS5cuYeHChfjqq6/QvXt3BAUFYebMmYiOjtbUycrKwsWLF9GvXz+djrN27Vrk5OTAx8cHP//8M27fvo0TJ07gxRdfRNOmTbFw4cIyt584cSJu3LiBgwcP6nS8p8FuGqq1FGam6DFrnNRhEFEl8W7ojT4N+lT5DKxBQUGYNGkSXnnlFVhZWeHdd9/F7du3SxwHUpKXX34ZxsbGCA4Oho+PDx49egRfX1/4+/tj4MCBAIAJEybg6NGjGDlyJH7++WfI5XJ89913cHNzQ8+ePXU6TpMmTfDbb7/hgw8+wNChQ5GYmAghBIYMGYLt27fD3Ny8zO3t7OwwatQofPDBBxgyZAiMjCqv/UImqmp6taeQnp4Oa2trpKWlFRtZTERENc+jR48QExMDDw8Pnf+IV1eZmZmoV68eVq5cibFjx+q0zbp163D48GGcPHlS5+N069YNM2bMwPDhwysaKgIDAxEUFITTp0+jW7duFd7Pv5X1Wur695stI0RERHq4dOkSIiMjoVKpkJaWhsWLFwMABg8erPM+Jk6ciNTUVGRkZOh8b5ohQ4bgzTffrHDcALBo0SK4u7vjl19+gUqlqtTWDn0wGSEiItLTihUrEBUVBYVCgU6dOiE0NBT29vY6b29sbFxsBtSy2Nvb491339Usx8XFoWXLlqXWv3btGtzc3EpcN3r0aJ2PW1WYjBAREemhQ4cOuHjxoqQxuLq6IiIiosz1NQmTESIiohrG2Ni41HlMaqLq0VlEREREzywmI0RERCQpdtMQGVh+QX6Vz3tARFSTMRkhMqDg2GAsC1umdXtzJ3MnzFXNrbQZIYmIajp20xAZSHBsMALOBCApIwFeNzzRP6IDvG544kFGIgLOBCA4NljqEImIqiW2jBAZQH5BPpaFLcOLf7RGs+TXkKuwBQC4ZQOd76cgyu4glpsvR58GfdhlQ2RAIj8fWb9fRF5SEowdHGDeuRNk8przGRs5ciRatGiB9957T6f6c+fORWZmJtasWVPJkVUttowQGUB4Yjja/eoIz4yxyDWx0VqXa2IDz4yxaPurA8ITw6UJkKgWSj91Cn/19Uacnx/uzZ6NOD8//NXXG+mnTkkWU5s2bTBp0qQS123fvh1KpRIPHjwAAFy+fBnHjh3DjBkzdN7/7NmzsW3bNq0b6JXl1q1bkMlkJc5JcubMGchkMqSmphZb5+7ujtWrV+sc19NiMkJkAPEp99As+bXChSfv3Pl4uVnyEMSn3KviyIhqp/RTp3B35lvIi4/XKs9LSMDdmW9JlpCMHTsWe/bsQXZ2drF1W7ZswaBBgzQzta5Zswavv/466tSpo/P+7e3t4ePjg88//9xgMVcHTEaIDCDn++uFXTOl3UJcJkOuwg4531+v2sCIaiGRn4+EJUuBku7z+rgsYclSiPz8Sjl+RkYGRowYAQsLC7i4uGDVqlXo3bs33nrrLfj6+iI7OxsHDx7U2iYmJgZnzpzR3EgvPz8fBw4c0NylFwAiIyNhbm6OXbt2acr27dsHMzMzXLt2TVM2cOBA7Nmzp1Kem1SYjBAZgFWybn3UutYjotJl/X6xWIuIFiGQFx+PrN8rZ8r2gIAAnDt3DocPH8bp06cRGhqK8PDCLlh7e3sMHjwYmzdv1tpm69atqF+/Pvr16wcAuHLlCtLS0tC5c2dNnebNm2PFihWYMmUK4uLicOfOHUyaNAnLly/Xug+NSqXCnTt3cOvWrUp5flLgAFYiAzB3sgVSdaxHRE8lLynJoPX0kZGRgW3btmHXrl3o27cvgMLul3/fC2bs2LF4+eWXERMTAw8PDwghsG3bNvj5+WnukhsbGwu5XA5HR0et/U+ZMgXHjh2Dr68vFAoFunTpgunTp2vVKTpWbGws3N3dDf4cpcCWESIDUE3yhTInpeRmYwAQAsqcZKgm+VZtYES1kLGDg0Hr6SM6Ohq5ublQqVSaMmtrazRr1kyz/OKLL6J+/frYsmULACAkJARxcXFad8vNzs6GUqmErISu3c2bN+PKlSsIDw/H1q1bi9UxMzMDAGRlZRn0uUmJyQiRASjMTOHhFle48GRC8njZw+02FGamVRwZUe1j3rkTjJ2dyxyjZezsDPPOnao2sMeMjIzg7++Pbdu2oaCgAFu2bEGfPn3g6empqWNvb4+srCyo1epi21++fBmZmZnIzMzE/fv3i61PTk4GADg8ZbJlZWUFAEhLSyu2LjU1FdbW1k+1f31UKBlZt24d3N3dYWpqiq5duyIsLKzUurm5uVi8eDEaNWoEU1NTtGvXDidOnKhwwETVVd9Fs9Dc+QqU6lStcqU6Bc2dr6DvolnSBEZUy8jkcji9N+/xQslXrzm9N69S5hvx9PSEiYkJfvvtN01ZWloabty4oVVv9OjRuH37Ng4dOoRvvvlGM3C1SPv27QFAa2AqUJho+Pv7Y/78+fD398eIESOKXZlz9epVmJiYoFWrVk/1XJo0aQIjIyNcvKg9tiY6OhppaWlo2rTpU+1fH3qPGdm7dy8CAgKwYcMGdO3aFatXr4aPjw+ioqKK9X0BwIIFC7Bjxw5s3LgRzZs3x8mTJ/Gf//wH58+fR4cOHQzyJIiqi76LZkGd/QhhG3YgKzEV5o42UE3yZYsIkYFZ9esHfLoaCUuWag1mNXZygtN78wrXVwJLS0v4+fnhnXfegZ2dHRwdHREYGAgjIyOt7hQPDw+88MILmDBhApRKJYYMGaK1HwcHB3Ts2BFnz57VJCYAMGnSJDRo0AALFixATk4OOnTogNmzZ2PdunWaOqGhoejZs6emu0YXUVFRxcpatWqFcePG4e2334axsTHatGmD27dvY86cOejWrRuee+45Pc7MUxJ6UqlUYurUqZrl/Px84erqKpYuXVpifRcXF7F27VqtsiFDhogRI0bofMy0tDQBQKSlpekbLhERVUPZ2dni2rVrIjs7+6n2U5CXJ/755VeR+v0R8c8vv4qCvDwDRVi69PR0MXz4cGFubi6cnZ1FUFCQUKlUYu7cuVr1du3aJQCIKVOmlLif9evXi27dummWt23bJiwsLMSNGzc0Zb/++qswMTERx44d05Q1a9ZM7N69W6dYY2JiBIASH7dv3xbZ2dkiMDBQNG/eXJiZmQkPDw8xYcIEkZSUpPP5KOu11PXvt0yI0kbcFadWq2Fubo4DBw7g1Vdf1ZT7+fkhNTUV3333XbFt6tati48//liricrX1xdnz54t9bKknJwc5OTkaJbT09PRoEEDpKWlafq4iIio5nr06JHmahNT05rdcpiZmYl69eph5cqVxbpjypKdnY1mzZph79698PLy0mmb48eP4+2338aVK1dgbFw9Logt67VMT0+HtbV1uX+/9Roz8uDBA+Tn58PJyUmr3MnJCfGlXPPt4+ODoKAg3Lx5EwUFBTh9+jQOHTpU4qCcIkuXLoW1tbXm0aBBA33CJCIiqjSXLl3C7t278ffffyM8PBwjRowAAAwePFiv/ZiZmeHrr7/WTA+vi8zMTGzZsqXaJCKGUulX03z66ado0qQJmjdvDoVCgWnTpmH06NGaa61LMm/ePKSlpWket2/fruwwiYiIdLZixQq0a9cO3t7eyMzMRGhoqGaad3307t1baxbW8vz3v/9F165dNcuTJk1CnTp1SnyUdo+c6kiv1Mre3h5yuRwJCQla5QkJCXB2di5xGwcHB3z77bd49OgRHj58CFdXV8ydO1frEqcnKZVKKJVKfUIjIiKqEh06dCh2BYpUFi9ejNmzZ5e4riYNa9ArGVEoFOjUqRNCQkI0Y0YKCgoQEhKCadOmlbmtqakp6tWrh9zcXBw8eBBDhw6tcNBEREQEODo6lngla02jd6dTQEAA/Pz80LlzZ6hUKqxevRqZmZmameVGjRqFevXqYenSpQCAX3/9FXfv3kX79u1x9+5dfPDBBygoKMC7775r2GdCRERENZLeyciwYcOQlJSEhQsXIj4+Hu3bt8eJEyc0g1rj4uK0xoM8evQICxYsQHR0NOrUqYP+/ftj+/btsLGxMdiTICIioppLr0t7paLrpUFERFQz1KZLe591VX5pLxEREZGhMRkhIiIiSTEZISIiqmJRUVFwdnZGRkaGTvUfPHgAR0dH3Llzp5IjkwaTESIiqrEKCgTuRqXgxm/xuBuVgoIC6YdB9u7dGzKZrNgjLy9PU2fevHmYPn06LC0tddqnvb09Ro0ahcDAQJ3j8Pf317p1y5MxvvXWW8XKt27dKskFJrVrPlmiZ1h+QT7CE8ORlJUEB3MHdHTsCLmR4W+hTlRd/H0pEaF7byIz9f/vZWZho0TPYU3QqIO0c2+MHz8eixcv1iormsI9Li4OR44cwZo1a/Ta5+jRo9GpUyd88sknsLOzM1is1QFbRohqgeDYYLy0zwcb1/wPP63Zg41r/oeX9vkgODZY6tCIKsXflxJx4ourWokIAGSm5uDEF1fx96XESjt2RkYGRowYAQsLC7i4uGDVqlXFWhrMzc3h7Oys9Siyb98+tGvXDvXq1dOUjRkzBm3bttXcJFatVqNDhw4YNWqUpk6rVq3g6uqKb775ptKem1SYjBDVcMGxwTi+ZjNePzUZ7R7OhFu2P9o9nInXT03G8TWbmZBQrVNQIBC692aZdc7uu1lpXTYBAQE4d+4cDh8+jNOnTyM0NBTh4eE6bx8aGorOnTtrlX322WfIzMzE3LlzAQDz589Hamoq1q5dq1VPpVIhNDT06Z9ENcNkhKgGyy/Ix49rtsMzYyxyTWy01uWa2MAzYyx+XLMD+QX50gRIVAnu30wt1iLypH9ScnD/ZqrBj52RkYFt27ZhxYoV6Nu3L1q3bo0tW7YgP1/7M7Z+/Xqtm9a9/fbbmnWxsbFwdXXVql+nTh3s2LED69atw8KFC7F69Wps37692Nwcrq6uiI2NNfjzkhrHjBDVYGFxv6Bh8iDkmgCQybRXymSAEGiY/ArC4n6Bl3t3SWIkMrTM9LITEX3r6SM6Ohq5ublQqVSaMmtrazRr1kyr3ogRIzB//nzN8r8HhWZnZ5c40ZuXlxdmz56NDz/8EHPmzEGPHj2K1TEzM0NWVpYBnkn1wpYRohrsrz3ByFXYFk9EishkyFXY4a897Kqh2sPCSre7uutarzJYW1ujcePGmoe9vb1mnb29PVJSUoptU1BQgHPnzkEul+Ovv/4qcb/JyclwcHB46visrKyQlpZWrDw1NRXW1tZPvX99MRkhqsEUybr1ietaj6gmcGliAwubshONOrZKuDSxMfixPT09YWJigt9++01TlpaWhhs3bui8jw4dOuDatWvFyj/55BNERkbip59+wokTJ7Bly5Zida5evYoOHTpULPh/adasWYnjXMLDw9G0adOn3r++mIwQ1WD1nOuVX0mPekQ1gZGRDD2HNSmzTo+hTWBkVEqL4VOwtLSEn58f3nnnHfz444/4888/MXbsWBgZGUFWWgvlE3x8fHDhwgWtcSaXLl3CwoUL8dVXX6F79+4ICgrCzJkzER0dramTlZWFixcvol+/fjrHm5aWhoiICK3H7du3MXnyZNy4cQMzZszAlStXEBUVhaCgIOzevVtrfEtVYTJCVIOpJo+EIicFKO1+l0JAkZMM1eSRVRsYUSVr1MERL01sXayFpI6tEi9NbF2p84wEBQXBy8sLr7zyCry9vdG9e3e0aNFC5xv+vfzyyzA2NkZwcGH36aNHj+Dr6wt/f38MHDgQADBhwgT06dMHI0eO1CQt3333Hdzc3NCzZ0+dYz1z5gw6dOig9Vi0aBE8PT3x888/IzIyEt7e3ujatSv27duH/fv346WXXtLzjDw93rWXqIYLCVyFyPi2hQv//mX2+KPd3PkK+i6aJUFkRKUz1F17CwpE4dU16TmwsCrsmqmMFpGyZGZmol69eli5ciXGjh2r0zbr1q3D4cOHcfLkSZ2P061bN8yYMQPDhw+vaKiVwhB37eXVNEQ1XN9Fs4DAVYiJc0OO0lZTrlSnwMPtNhMRqtWMjGSo18y2/IoGdOnSJURGRkKlUiEtLU0z0+rgwYN13sfEiRORmpqKjIwMnaaEf/DgAYYMGYI333yzwnFXZ0xGiGqBvotmQZ39CGEbdiArMRXmjjZQTfKFwqzivziJqHQrVqxAVFQUFAoFOnXqhNDQUK0rZspjbGysdelveezt7fHuu+9qluPi4tCyZctS61+7dg1ubm46719qTEaIagmFmSl6zBondRhEtV6HDh1w8eJFSWNwdXVFREREmetrEiYjRERENYyxsTEaN24sdRgGw6tpiIhIMjXgGgoqhyFeQyYjRERU5UxMTACgVk5t/qwpeg2LXtOKYDcNERFVOblcDhsbGyQmJgIAzM3NdZ40jKoHIQSysrKQmJgIGxsbyOXyCu+LyQgREUnC2dkZADQJCdVMNjY2mteyopiMEBGRJGQyGVxcXODo6Ijc3Fypw6EKMDExeaoWkSJMRoiISFJyudwgf9Co5uIAViIiIpIUkxEiIiKSFJMRIiIikhSTESIiIpIUkxEiIiKSFJMRIiIikhSTESIiIpIUkxEiIiKSFCc9I6Knos5+hLANO5CVmApzRxuoJvlCYWYqdVhEVIMwGSGiCgsJXIWYODfkKD0LC1KAyClH4eEWh76LZkkbHBHVGOymIaIKCQlchcj4tshR2GiV5yhsEBnfFiGBq6QJjIhqHCYjRKQ3dfYjRMe5FS48edv3x8vRcQ2gzn5UxZERUU3EZISI9Bb2+XaolbbFE5EiMhnUSjuEfb69agMjohqJyQgR6e1u/F2D1iOiZxuTESLSm9qulBaRCtYjomcbkxEi0lvjN7xhok4BhCi5ghAwUSej8RveVRsYEdVITEaISG8qt26ItTtcuPBkQvJ4OdbuCFRu3ao4MiKqiZiMEJHe5EZy9Jk+EtGWm2CSm6q1ziQ3BdGWm9Bnui/kRnJpAiSiGoWTnhFRhXg39AamA8svLIPHHxawzbJGinkaYtpkYo7X3ML1REQ6kAlRWqdv9ZGeng5ra2ukpaXByspK6nCI6F/yC/IRnhiOpKwkOJg7oKNjR7aIEBEA3f9+s2WEiJ6K3EiOLs5dpA6DiGowjhkhIiIiSTEZISIiIkkxGSEiIiJJMRkhIiIiSTEZISIiIkkxGSEiIiJJMRkhIiIiSTEZISIiIkkxGSEiIiJJMRkhIiIiSTEZISIiIkkxGSEiIiJJMRkhIiIiSTEZISIiIkkxGSEiIiJJGUsdABFRRaizHyFsww5kJabC3NEGqkm+UJiZSh0WEVVAhVpG1q1bB3d3d5iamqJr164ICwsrs/7q1avRrFkzmJmZoUGDBpg1axYePXpUoYCJiEICV+HrKUdxOcoTN1M64nKUJ76echQhgaukDo2IKkDvZGTv3r0ICAhAYGAgwsPD0a5dO/j4+CAxMbHE+rt27cLcuXMRGBiI69evY9OmTdi7dy/ee++9pw6eiJ49IYGrEBnfFjkKG63yHIUNIuPbMiEhqoH0TkaCgoIwfvx4jB49Gi1btsSGDRtgbm6OzZs3l1j//Pnz6N69O4YPHw53d3f069cPb775ZrmtKURET1JnP0J0nFvhgkymvfLxcnRcA6iz2fJKVJPolYyo1WpcvHgR3t7e/78DIyN4e3vjwoULJW7z3HPP4eLFi5rkIzo6GseOHUP//v1LPU5OTg7S09O1HkREYZ9vh1ppWzwRKSKTQa20Q9jn26s2MCJ6KnoNYH3w4AHy8/Ph5OSkVe7k5ITIyMgStxk+fDgePHiAHj16QAiBvLw8TJo0qcxumqVLl2LRokX6hEZEz4C78XcBNNKxHhHVFJV+ae+ZM2ewZMkSrF+/HuHh4Th06BCOHj2KDz/8sNRt5s2bh7S0NM3j9u3blR0mEdUAartSWkQqWI+Iqge9Wkbs7e0hl8uRkJCgVZ6QkABnZ+cSt3n//fcxcuRIjBs3DgDQpk0bZGZmYsKECZg/fz6MjIrnQ0qlEkqlUp/QiOgZ0PgNb/yxOB65JjYld9UIAZPcFDR+w7v4OiKqtvRqGVEoFOjUqRNCQkI0ZQUFBQgJCYGXl1eJ22RlZRVLOORyOQBACKFvvET0DFO5dUOs3eHChSe/Px4vx9odgcqtWxVHRkRPQ+9umoCAAGzcuBHbtm3D9evXMXnyZGRmZmL06NEAgFGjRmHevHma+gMHDsTnn3+OPXv2ICYmBqdPn8b777+PgQMHapISIiJdyI3k6DN9JKItN8EkN1VrnUluCqItN6HPdF/IjfjdQlST6D0D67Bhw5CUlISFCxciPj4e7du3x4kTJzSDWuPi4rRaQhYsWACZTIYFCxbg7t27cHBwwMCBA/HRRx8Z7lkQ0TPDu6E3MB1YfmEZPP6wgG2WNVLM0xDTJhNzvOYWrieiGkUmakBfSXp6OqytrZGWlgYrKyupwyGiaiC/IB/hieFIykqCg7kDOjp2ZIsIUTWj699v3puGiGokuZEcXZy7SB0GERkA79pLREREkmIyQkRERJJiMkJERESSYjJCREREkmIyQkRERJJiMkJERESSYjJCREREkmIyQkRERJJiMkJERESSYjJCREREkmIyQkRERJJiMkJERESSYjJCREREkmIyQkRERJJiMkJERESSMpY6ACKimiS/IB/hieFIykqCg7kDOjp2hNxILnVYRDUakxEiIh0FxwZj+YVl8PjDArZZ1kgxT0NMm0zM8ZoL74beUodHVGMxGSEi0kFwbDCOr9mM15MnI1dhCwBwywY6n0rB8d83A9PBhISogjhmhIioHPkF+fhxzXZ4ZoxFromN1rpcExt4ZozFj2t2IL8gX5oAiWo4JiNEROUIi/sFDZMHFS7IZNorHy83TH4FYXG/VHFkRLUDkxEionL8tSe4sGvmyUSkiEyGXIUd/toTXLWBEdUSTEaIiMqhSBYGrUdE2piMEBGVo55zPYPWIyJtTEaIiMqhmjwSipwUQJTS8iEEFDnJUE0eWbWBEdUSTEaIiMqhMDOFp1tc4cKTCcnjZU+321CYmVZxZES1A5MRIiId9F00C82dr0CpTtUqV6pT0Nz5CvoumiVNYES1ACc9IyLSUd9Fs6DOfoSwDTuQlZgKc0cbqCb5skWE6CkxGSEi0oPCzBQ9Zo2TOgyiWoXdNERERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCkmI0RERCQpJiNEREQkKSYjREREJCljqQMgIqLy5RfkIzwxHElZSXAwd0BHx46QG8mlDovIIJiMEBFVc8GxwVgWtgwJWQmaMidzJ8xVzYV3Q28JIyMyDHbTEBFVY8GxwQg4E4CkjAR43fBE/4gO8LrhiQcZiQg4E4Dg2GCpQyR6amwZISKqpvIL8rEsbBle/KM1miW/hlyFLQDALRvofD8FUXYHsdx8Ofo06MMuG6rR2DJCRFRNhSeGo92vjvDMGItcExutdbkmNvDMGIu2vzogPDFcmgCJDITJCBFRNRWfcg/Nkl8rXJDJtFc+Xm6WPATxKfeqODIiw2IyQkRUTeV8f72wa+bJRKSITIZchR1yvr9etYERGRiTESKiasoqWbdxILrWI6qumIwQEVVT5k62Bq1HVF0xGSEiqqZUk3yhzEkBhCi5ghBQ5iRDNcm3agMjMjAmI0RE1ZTCzBQebnGFC08mJI+XPdxuQ2FmWsWRERkWkxEiomqs76JZaO58BUp1qla5Up2C5s5X0HfRLGkCIzIgTnpGRFTN9V00C+rsRwjbsANZiakwd7SBapIvW0So1mAyQkRUAyjMTNFj1jipwyCqFOymISIiIkkxGSEiIiJJVSgZWbduHdzd3WFqaoquXbsiLCys1Lq9e/eGTCYr9hgwYECFgyYiIqLaQ+9kZO/evQgICEBgYCDCw8PRrl07+Pj4IDExscT6hw4dwv379zWPq1evQi6X4/XXX3/q4ImIiKjm0zsZCQoKwvjx4zF69Gi0bNkSGzZsgLm5OTZv3lxifTs7Ozg7O2sep0+fhrm5OZMRIiIiAqBnMqJWq3Hx4kV4e3v//w6MjODt7Y0LFy7otI9NmzbhjTfegIWFRal1cnJykJ6ervUgIiKi2kmvZOTBgwfIz8+Hk5OTVrmTkxPi4+PL3T4sLAxXr17FuHFlX562dOlSWFtbax4NGjTQJ0wiIiKqQar0appNmzahTZs2UKlUZdabN28e0tLSNI/bt29XUYRERERU1fSa9Mze3h5yuRwJCQla5QkJCXB2di5z28zMTOzZsweLFy8u9zhKpRJKpVKf0IiIiKiG0qtlRKFQoFOnTggJCdGUFRQUICQkBF5eXmVuu3//fuTk5MDXl3eXJCIiov+n93TwAQEB8PPzQ+fOnaFSqbB69WpkZmZi9OjRAIBRo0ahXr16WLp0qdZ2mzZtwquvvoq6desaJnIiIiKqFfRORoYNG4akpCQsXLgQ8fHxaN++PU6cOKEZ1BoXFwcjI+0Gl6ioKJw9exanTp0yTNRERERUa8iEEELqIMqTnp4Oa2trpKWlwcrKSupwiIiISAe6/v3mvWmIiIhIUkxGiIiISFJMRoiIiEhSeg9gJSKiZ0d2ZiaOrluD7Af/wMy+DgZMnQ6zMm7nQVQRTEaIiKhEO+cvROb9dshVFM6anZ4ObJ9+AhYulzHio/InsCTSFbtpiIiomJ3zFyL1QS/kmtholeea2CD1QS/snL9QmsCoVmIyQkREWrIzM5F5v13hgkymvfLxcub9tsjOzKziyKi2YjJCRERajq5bg1yFbfFEpIhMhlyFHY6uW1O1gVGtxWSEiIi0ZD/4x6D1iMrDZISIiLSY2dcxaD2i8jAZISIiLQOmToeJOgUo7W4hQsBEnYwBU6dXbWBUazEZISIiLWYWFrBwuVy48GRC8njZwuUK5xshg2EyQkRExYz4aDFs7H+CSW6qVrlJbgps7H/iPCNkULxrLxERlYozsNLT0PXvN5MRIiIiqhS6/v1mNw0RERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJyljqAIiIiJ6kzn6EsA07kJWYCnNHG6gm+UJhZip1WFRJmIwQEVG1EhK4CjFxbshRehYWpACRU47Cwy0OfRfNkjY4qhTspiEiomojJHAVIuPbIkdho1Weo7BBZHxbhASukiYwqlRMRoiIqFpQZz9CTJxb4YJMpr3y8XJMXAOosx9VcWRU2ZiMEBFRtRC2YQdylLbFE5EiMhlylHYI27CjagOjSsdkhIiIqoWshBSD1qOag8kIERFVC+l2+QatRzUHkxEiIqoWlANbwESdAghRcgUhYKJOhnJgi6oNjCodkxEiIqoWnG1dEWV3sHDhyYTk8XKU3SE427pWcWRU2ZiMEBFRtdDRsSMud01EtOUmmOSmaq0zyU1BtOUmXOmahI6OHaUJkCoNJz0jIqJqQW4kx1zVXARkBSAk/ypUf3vANssaKeZpCGsUgwI5EKQKgtxILnWoZGAyIUrrnKs+0tPTYW1tjbS0NFhZWUkdDhERVaLg2GAsC1uGhKwETZmzuTPmqObAu6G3hJGRvnT9+82WESIiqla8G3qjT4M+CE8MR1JWEhzMHdDRsSNbRGoxJiNERFTtyI3k6OLcReowqIpwACsRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJiskIERERSYrJCBEREUmKyQgRERFJinftJSIiekyd/QhhG3YgKzEV5o42UE3yhcLMVOqwaj0mI0RERABCAlchJs4NOUrPwoIUIHLKUXi4xaHvolnSBlfLsZuGiIieeSGBqxAZ3xY5Chut8hyFDSLj2yIkcJU0gT0jKpSMrFu3Du7u7jA1NUXXrl0RFhZWZv3U1FRMnToVLi4uUCqVaNq0KY4dO1ahgImIiAxJnf0I0XFuhQsymfbKx8vRcQ2gzn5UxZE9O/RORvbu3YuAgAAEBgYiPDwc7dq1g4+PDxITE0usr1ar8eKLL+LWrVs4cOAAoqKisHHjRtSrV++pgyciInpaYZ9vh1ppWzwRKSKTQa20Q9jn26s2sGeI3mNGgoKCMH78eIwePRoAsGHDBhw9ehSbN2/G3Llzi9XfvHkzkpOTcf78eZiYmAAA3N3dny5qIiIiA7kbfxdAIx3rUWXQq2VErVbj4sWL8Pb2/v8dGBnB29sbFy5cKHGbw4cPw8vLC1OnToWTkxNat26NJUuWID8/v9Tj5OTkID09XetBRERUGdR2pbSIVLAe6U+vZOTBgwfIz8+Hk5OTVrmTkxPi4+NL3CY6OhoHDhxAfn4+jh07hvfffx8rV67E//73v1KPs3TpUlhbW2seDRo00CdMIiIinTV+wxsm6hRAiJIrCAETdTIav+Fd8np6apV+NU1BQQEcHR3x5ZdfolOnThg2bBjmz5+PDRs2lLrNvHnzkJaWpnncvn27ssMkIqJnlMqtG2LtDhcuPJmQPF6OtTsClVu3Ko7s2aFXMmJvbw+5XI6EhASt8oSEBDg7O5e4jYuLC5o2bQq5XK4pa9GiBeLj46FWq0vcRqlUwsrKSutBRERUGeRGcvSZPhLRlptgkpuqtc4kNwXRlpvQZ7ov5EbykndAT02vAawKhQKdOnVCSEgIXn31VQCFLR8hISGYNm1aidt0794du3btQkFBAYyMCnOfGzduwMXFBQqF4umiJyIiMgDvht7AdGD5hWXw+MMCtlnWSDFPQ0ybTMzxmlu4niqN3lfTBAQEwM/PD507d4ZKpcLq1auRmZmpubpm1KhRqFevHpYuXQoAmDx5MtauXYuZM2di+vTpuHnzJpYsWYIZM2YY9pkQERE9Be+G3ujToA/Ce4cjKSsJDuYO6OjYkS0iVUDvZGTYsGFISkrCwoULER8fj/bt2+PEiROaQa1xcXGaFhAAaNCgAU6ePIlZs2ahbdu2qFevHmbOnIk5c+YY7lkQEREZgNxIji7OXaQO45kjE6K04cPVR3p6OqytrZGWlsbxI0RERDWErn+/eW8aIiIikhSTESIiIpIUkxEiIiKSFJMRIiIikhSTESIiIpIUkxEiIiKSFJMRIiIikhSTESIiIpIUkxEiIiKSlN7TwRMREZFh5RfkIzzx2b0nDpMRIiIiCQXHBmNZ2DIkZCVoypzMnTBX9ezcLZjdNERERBIJjg1GwJkArUQEABKzEhFwJgDBscESRVa1mIwQERFJIL8gH8vClkFAwCgf8Lrhif4RHeB1wxOy/MJ72C4PW478gnyJI6187KYhIiKSQHhiOBKyEtDvjzZolvwachW2AAC3bKDz/RRE2R3EqTZ/IDwxHF2cu0gcbeViywgREZEEkrKS0O+PNvDMGItcExutdbkmNvDMGIt+f7RBUlaSNAFWISYjREREErA1skaz5NcKF2Qy7ZWPl5slD4GtkXUVR1b1mIwQERFJIP/QtcKumScTkSIyGXIVdsg/dK1qA5MAkxEiIiIJZCWlGbReTcZkhIiISALmjjYGrVeTMRkhIiKSgGqSL5Q5KYAQJVcQAsqcZKgm+VZtYBJgMkJERCQBhZkpPNziCheeTEgeL3u43YbCzLSKI6t6TEaIiIgk0nfRLDR3vgKlOlWrXKlOQXPnK+i7aJY0gVUxTnpGREQkob6LZkGd/QhhG3YgKzEV5o42UE3yfSZaRIowGSEiIpKYwswUPWaNkzoMybCbhoiIiCTFZISIiIgkxWSEiIiIJMVkhIiIiCTFZISIiIgkxWSEiIiIJMVkhIiIiCTFZISIiIgkxWSEiIiIJMVkhIiIiCTFZISIiIgkxWSEiIiIJMUb5RERET2jqsvdgpmMEBERPYNCAlchJs4NOUrPwoIUIHLKUXi4xaHvollVGgu7aYiIiJ4xIYGrEBnfFjkKG63yHIUNIuPbIiRwVZXGw2SEiIjoGaLOfoToOLfCBZlMe+Xj5ei4BlBnP6qymJiMEBERPUPCPt8OtdK2eCJSRCaDWmmHsM+3V1lMTEaIiIieIXfj7xq0niEwGSEiInqGqO1KaRGpYD1DYDJCRET0DGn8hjdM1CmAECVXEAIm6mQ0fsO7ymJiMkJERPQMUbl1Q6zd4cKFJxOSx8uxdkegcutWZTExGSEiInqGyI3k6DN9JKItN8EkN1VrnUluCqItN6HPdF/IjeRVFhMnPSMiInrGeDf0BqYDyy8sg8cfFrDNskaKeRpi2mRijtfcwvVVSCZEaZ1G1Ud6ejqsra2RlpYGKysrqcMhIiKqFfIL8hGeGI6krCQ4mDugo2NHg7aI6Pr3my0jREREzyi5kRxdnLtIHQbHjBAREZG0mIwQERGRpJiMEBERkaSYjBAREZGkmIwQERGRpJiMEBERkaSYjBAREZGkmIwQERGRpJiMEBERkaRqxAysRTPWp6enSxwJERER6aro73Z5d56pEclIRkYGAKBBgwYSR0JERET6ysjIgLW1danra8SN8goKCnDv3j1YWlpCJpMZbL/p6elo0KABbt++zRvwSYDnX3p8DaTF8y8tnv/KJ4RARkYGXF1dYWRU+siQGtEyYmRkhPr161fa/q2srPhGlBDPv/T4GkiL519aPP+Vq6wWkSIcwEpERESSYjJCREREknqmkxGlUonAwEAolUqpQ3km8fxLj6+BtHj+pcXzX33UiAGsREREVHs90y0jREREJD0mI0RERCQpJiNEREQkKSYjREREJCkmI0RERCSpWp+MrFu3Du7u7jA1NUXXrl0RFhZWZv39+/ejefPmMDU1RZs2bXDs2LEqirR20uf8b9y4ET179oStrS1sbW3h7e1d7utF5dP3M1Bkz549kMlkePXVVys3wFpO3/OfmpqKqVOnwsXFBUqlEk2bNuX30FPQ9/yvXr0azZo1g5mZGRo0aIBZs2bh0aNHVRTtM0zUYnv27BEKhUJs3rxZ/Pnnn2L8+PHCxsZGJCQklFj/3LlzQi6Xi48//lhcu3ZNLFiwQJiYmIg//vijiiOvHfQ9/8OHDxfr1q0Tly5dEtevXxf+/v7C2tpa3Llzp4ojrz30fQ2KxMTEiHr16omePXuKwYMHV02wtZC+5z8nJ0d07txZ9O/fX5w9e1bExMSIM2fOiIiIiCqOvHbQ9/zv3LlTKJVKsXPnThETEyNOnjwpXFxcxKxZs6o48mdPrU5GVCqVmDp1qmY5Pz9fuLq6iqVLl5ZYf+jQoWLAgAFaZV27dhUTJ06s1DhrK33P/5Py8vKEpaWl2LZtW2WFWOtV5DXIy8sTzz33nPjqq6+En58fk5GnoO/5//zzz4Wnp6dQq9VVFWKtpu/5nzp1qnjhhRe0ygICAkT37t0rNU4SotZ206jValy8eBHe3t6aMiMjI3h7e+PChQslbnPhwgWt+gDg4+NTan0qXUXO/5OysrKQm5sLOzu7ygqzVqvoa7B48WI4Ojpi7NixVRFmrVWR83/48GF4eXlh6tSpcHJyQuvWrbFkyRLk5+dXVdi1RkXO/3PPPYeLFy9qunKio6Nx7Ngx9O/fv0pifpbViLv2VsSDBw+Qn58PJycnrXInJydERkaWuE18fHyJ9ePj4ystztqqIuf/SXPmzIGrq2uxBJF0U5HX4OzZs9i0aRMiIiKqIMLarSLnPzo6Gj/88ANGjBiBY8eO4a+//sKUKVOQm5uLwMDAqgi71qjI+R8+fDgePHiAHj16QAiBvLw8TJo0Ce+9915VhPxMq7UtI1SzLVu2DHv27ME333wDU1NTqcN5JmRkZGDkyJHYuHEj7O3tpQ7nmVRQUABHR0d8+eWX6NSpE4YNG4b58+djw4YNUof2TDhz5gyWLFmC9evXIzw8HIcOHcLRo0fx4YcfSh1arVdrW0bs7e0hl8uRkJCgVZ6QkABnZ+cSt3F2dtarPpWuIue/yIoVK7Bs2TIEBwejbdu2lRlmrabva/D333/j1q1bGDhwoKasoKAAAGBsbIyoqCg0atSocoOuRSryGXBxcYGJiQnkcrmmrEWLFoiPj4darYZCoajUmGuTipz/999/HyNHjsS4ceMAAG3atEFmZiYmTJiA+fPnw8iIv98rS609swqFAp06dUJISIimrKCgACEhIfDy8ipxGy8vL636AHD69OlS61PpKnL+AeDjjz/Ghx9+iBMnTqBz585VEWqtpe9r0Lx5c/zxxx+IiIjQPAYNGoQ+ffogIiICDRo0qMrwa7yKfAa6d++Ov/76S5MEAsCNGzfg4uLCRERPFTn/WVlZxRKOosRQ8J6ylUvqEbSVac+ePUKpVIqtW7eKa9euiQkTJggbGxsRHx8vhBBi5MiRYu7cuZr6586dE8bGxmLFihXi+vXrIjAwkJf2PgV9z/+yZcuEQqEQBw4cEPfv39c8MjIypHoKNZ6+r8GTeDXN09H3/MfFxQlLS0sxbdo0ERUVJY4cOSIcHR3F//73P6meQo2m7/kPDAwUlpaWYvfu3SI6OlqcOnVKNGrUSAwdOlSqp/DMqNXJiBBCrFmzRri5uQmFQiFUKpX45ZdfNOt69eol/Pz8tOrv27dPNG3aVCgUCtGqVStx9OjRKo64dtHn/Dds2FAAKPYIDAys+sBrEX0/A//GZOTp6Xv+z58/L7p27SqUSqXw9PQUH330kcjLy6viqGsPfc5/bm6u+OCDD0SjRo2EqampaNCggZgyZYpISUmp+sCfMTIh2PZERERE0qm1Y0aIiIioZmAyQkRERJJiMkJERESSYjJCREREkmIyQkRERJJiMkJERESSYjJCREREkmIyQkRERJJiMkJERESSYjJCREREkmIyQkRERJL6P6FxcdlFKd6pAAAAAElFTkSuQmCC"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the g_V and g_F when M = 16, N = 4, 8\n",
    "def plot_least_squares_qr(m, n):\n",
    "    x = np.linspace(0, 1, m + 1)[:-1]\n",
    "    y = 1 / (1 + x**2)\n",
    "    vm = vandermonde_matrix_mn(m, n)\n",
    "    fm = fourier_basis_matrix_mn(m, n)\n",
    "\n",
    "    # solve the least square system with QR decomposition by minimizing the 2-norm of the residual\n",
    "    vm_q, vm_r = np.linalg.qr(vm)\n",
    "    fm_q, fm_r = np.linalg.qr(fm)\n",
    "    vm_qt = np.transpose(vm_q)\n",
    "    fm_qt = np.transpose(fm_q)\n",
    "    vm_b = np.dot(vm_qt, y)\n",
    "    fm_b = np.dot(fm_qt, y)\n",
    "    vm_x = sp.linalg.solve_triangular(vm_r, vm_b)\n",
    "    fm_x = sp.linalg.solve_triangular(fm_r, fm_b)\n",
    "    vm_y = np.dot(vm, vm_x)\n",
    "    fm_y = np.dot(fm, fm_x)\n",
    "    plt.plot(x, y, 'o', label='f(x)')\n",
    "    plt.plot(x, vm_y, 'o', label='gV(x)_QR')\n",
    "    plt.plot(x, fm_y, 'o', label='gF(x)_QR')\n",
    "\n",
    "    # solve the least square system with LU decomposition by minimizing the 2-norm of the residual\n",
    "    x = np.linspace(0, 1, m + 1)[:-1]\n",
    "    y = 1 / (1 + x**2)\n",
    "    vm = vandermonde_matrix(m)\n",
    "    fm = fourier_basis_matrix(m)\n",
    "    vm_lu = sp.linalg.lu_factor(vm)\n",
    "    fm_lu = sp.linalg.lu_factor(fm)\n",
    "    vm_x = sp.linalg.lu_solve(vm_lu, y)\n",
    "    fm_x = sp.linalg.lu_solve(fm_lu, y)\n",
    "    vm_y = np.dot(vm, vm_x)\n",
    "    fm_y = np.dot(fm, fm_x)\n",
    "    plt.plot(x, vm_y, 'o', label='gV(x)_LU')\n",
    "    plt.plot(x, fm_y, 'o', label='gF(x)_LU')\n",
    "    plt.legend()\n",
    "    plt.title('Least Squares with QR and LU decomposition' + ' M = ' + str(m) + ' N = ' + str(n))\n",
    "    plt.show()\n",
    "plot_least_squares_qr(16, 4)\n",
    "plot_least_squares_qr(16, 8)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-05-27T12:24:55.192324700Z",
     "start_time": "2023-05-27T12:24:54.124882800Z"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
